In order to probe the phenomenon picked out by the first PC of our two PCA
analyses, we fit a series of models. First, we model F1 on the basis of
amplitude to see how changes in amplitude affect each vowel’s F1. This is
motivated by the thought that amplitude drives the common movement on F1
found by the first PC in our PCA analysis (developed in corpus_pca.Rmd).
Second, we investigate the sociolinguistic question of whether and to what extent changes in amplitude suggest that a speaker is coming to the end of a discrete topical unit of a monologue.
For both questions we will first apply simple linear and linear mixed models, before turning to more sophisticated GAMM models.
In order for this document to be independently understandable, we briefly run through the phenomenon of interest from the previous supplementary materials.
We first load the required libraries and define global variables.
# Tidyverse and friends
library(tidyverse)
library(broom)
library(glue)
library(patchwork)
# Animations
library(gganimate)
library(magick)
# Interactive plots
library(plotly)
# File management
library(here)
# Data scaling
library(scales)
# GAMMs
library(mgcv)
library(itsadug)
library(gratia)
# Linear Mixed Models
library(lme4)
library(optimx)
library(car)
# For variance inflation function
library(car)
# parallel computing - only used for `detectCores` function.
library(parallel)
# Global variables for plotting
vowel_colours_with_foot <- c(
START = "#00B0F6",
STRUT = "#F8766D",
LOT = "#00BF7D",
TRAP = "#FF62BC",
FOOT = "#966432",
KIT = "#39B600",
NURSE = "#00BFC4",
THOUGHT = "#E76BF3",
DRESS = "#9590FF",
FLEECE = "#D89000",
GOOSE = "#A3A500"
)
# Order = order at which these vowels appear on the right side of the main plot
# of the model of F1 by amplitude which is used in the paper. We don't reorder
# for each plot in these supplementaries.
vowels <- c(
"START", "STRUT", "LOT", "TRAP", "FOOT", "KIT", "NURSE", "THOUGHT",
"DRESS", "FLEECE", "GOOSE"
)
# Sometimes it is useful to split plots between high vowels and others.
# We define high and low vowels here for this purpose.
high_vowels <- c(
"DRESS",
"GOOSE",
"THOUGHT",
"FLEECE",
"NURSE"
)
front_vowels <- c(
"DRESS",
"FLEECE",
"NURSE",
"GOOSE",
"TRAP"
)
# Random seed set for reproducibility.
set.seed(5)
knitr::include_graphics(here('plots', 'PCA_with_amplitude_varplot.png'))
Figure 1.1: Variables plots from PCA analysis.
As depicted in Figure 1.1, PC1 of our PCA analysis for both 60 second and 240 second intervals reveals that F1s of each vowel move together with amplitude, and that this effect explains around 7.7% of the variance for the 60 second intervals and around 10.3% of the variance for the 240 second intervals.
In SM2_interval_representation.Rmd, examples of amplitude over the course of a
monologue were presented. These seemed to suggest that amplitude systematically
drops over time (e.g. the bottom panels of Figure
1.2).
knitr::include_graphics(here('plots', 'QB_NZ_F_369_combined.png'))
Figure 1.2: Amplitude over the course of monologue for 60 and 240 second itnervals.
We will be interested in whether changes in F1 over the course of a monologue are explained by changes in amplitude and whether there is some other effect which might explain systematic shifts in F1 over the course of a monologue.
Connection between SM4 and the paper: The models reported in the paper are developed in Section 2.3 and Section 3.2. The remainder of the document consists of assumption checks and exploration of alternative methods. The fact that alternative methods produce compatible results provides a ‘sanity check’ on the methods which we do report.
We load the filtered data.
qb_vowels <- read_rds(
here('processed_data', 'Quakebox_filtered.rds')
)
qb_vowels
We scale the variables so that they are comparable across speakers. Because we
are unable to control for different recording environments, we will scale
amplitude so we are dealing with relative amplitude within a monologue. We will
also scale speaker formant values. In both cases, we are simply z-scoring using
the base R scale function. We also scale articulation rate, pitch, and speaker
length both within and across speakers. Both forms of scaling will be useful for
fitting models. We also scale time so that each monologue has a value running
from 0 to 1, where 0 is the first token in the monologue and 1 is the last.
# scale time, collect speaker length.
qb_vowels <- qb_vowels %>%
group_by(Speaker) %>%
rename(
time = Target.segments.start,
art_rate = utterance.articulation.rate,
pitch = MeanPitch
) %>%
# Within speaker scaling.
mutate(
speaker_scaled_time = rescale(time, to = c(0, 1)),
speaker_length = max(time),
speaker_scaled_amp_max = scale(intensity_max),
speaker_scaled_art_rate = scale(art_rate),
speaker_scaled_pitch = scale(pitch)
) %>%
ungroup() %>%
# Across speaker scaling.
mutate(
scaled_art_rate = scale(art_rate),
scaled_pitch = scale(pitch),
scaled_length = scale(speaker_length)
# We don't scale amplitude across speakers as we can't control for recording
# variation.
)
# Scale formant data
qb_vowels <- qb_vowels %>%
group_by(Speaker, Vowel) %>%
mutate(
speaker_scaled_F1 = scale(F1_50),
speaker_scaled_F2 = scale(F2_50)
) %>%
# Remove rows with missing F1
filter(
!is.na(speaker_scaled_F1)
) %>%
ungroup() %>%
# We may use across-speaker scaled F1 as a response.
mutate(
scaled_F1 = scale(F1_50),
scaled_F1 = scale(F2_50)
)
We now look at the distributions of the key variables.
Our initial models probe the relationship between amplitude and F1. We will use both scaled and unscaled F1 values. We first look at the unscaled values (Figure 2.1). Inspection of the figure suggests that all are roughly normally distributed, although with quite different amounts of variance and moderate differences in their skewness.
qb_vowels %>%
ggplot(
aes(
x = F1_50
)
) +
geom_histogram(stat="density") +
facet_wrap(vars(Vowel)) +
labs(
title = "Unscaled F1 Distributions by Vowel"
)
Figure 2.1: Distributions of unscaled F1 values for each vowel.
Figure 2.2 shows the (predictable) consequences
of scaling these variables. All are now centred around zero and ranging between
-2.5 and 2.5. This is unsurprising because
we filtered out tokens with s.d. filtering
at -2.5 and 2.5 in preprocessing.Rmd.
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_F1
)
) +
geom_histogram( stat="density") +
facet_wrap(vars(Vowel)) +
labs(
title = "Scaled F1 Distributions by Vowel"
)
Figure 2.2: Distributions of scaled F1 values for each vowel.
Our main variable of interest is the maximum amplitude of the word in which the vowel token originates. Taking our amplitude values at the word level rather than the vowel token level helps us to generate reliable amplitude readings using Praat.1 Figure 2.3 shows the distribution of maximum amplitude.
(
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_amp_max
)
) +
geom_histogram(binwidth=0.1) +
labs(
title = "Scaled",
x = "Scaled max amplitude"
)
) +
(
qb_vowels %>%
ggplot(
aes(
x = intensity_max
)
) +
geom_histogram(binwidth=1)+
labs(
title = "Unscaled",
x = "Unscaled max amplitude"
)
) +
plot_annotation(
title = "Max word amplitude distributions"
)
Figure 2.3: Maximum amplitude, scaled and unscaled.
The automatically generated scale on the \(x\)-axis in Figure 2.3 suggests that we have a few outliers on the left tail. We look at these values by filtering for scaled maximum amplitude below \(-5\).
qb_vowels %>%
filter(
speaker_scaled_amp_max < -5
) %>%
select( # Remove non-informative variables.
-c("MatchId", "TargetId", "URL")
)
Interestingly, almost all of these have no pitch information. This suggests a tracking problem in Praat. We look at how common missing pitch information is in the dataframe.
qb_vowels %>%
filter(
is.na(pitch)
) %>%
select( # Remove non-informative variables.
-c("MatchId", "TargetId", "URL")
)
There are about 100 times more entries with missing pitch than there are entries with scaled amplitude < 0.5. We won’t delete all of these data points.
Instead, given that the arguments for adopting the 2.5 standard deviation cut off for formant values apply to amplitude values just as well. So we will apply the same cut off to the amplitude data.
qb_vowels <- qb_vowels %>%
filter(
abs(speaker_scaled_amp_max) <= 2.5
)
The resulting distribution for amplitude is depicted in Figure 2.4. No surprises here.
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_amp_max
)
) +
geom_histogram(binwidth=0.1) +
labs(
title = "Max amplitude distribution",
subtitle = "Filtered and scaled",
x = "Scaled max amplitude"
)
Figure 2.4: Filtered and scaled max amplitude distribution.
Our models will use articulation rate and pitch as control variables. When we use scaled F1 as the response variable, we will also scale these within speaker. We now check the distributions of scaled articulation rate and pitch.
qb_vowels %>%
ggplot(
aes(
x = scaled_pitch,
colour = participant_gender
)
) +
geom_freqpoly(binwidth=0.1, size=1, stat='density') +
labs(
title = "Across-speaker scaled mean pitch distribution",
x = "Scaled mean pitch"
)
Figure 2.5: Across speaker scaled mean pitch distribution.
Figure 2.5 shows that our across speaker scaling results in a bimodal distribution with modes mostly corresponding to the gender of the speaker. This scaling will only be used to aid the fitting process, so we don’t need to worry about it.
Figure 2.6 shows the distribution as a histogram of raw token counts.
qb_vowels %>%
ggplot(
aes(
x = scaled_pitch
)
) +
geom_histogram(binwidth=0.1) +
labs(
title = "Across-speaker scaled mean pitch distribution",
x = "Scaled mean pitch"
)
Figure 2.6: Histogram of across speaker scaled mean pitch distribution.
Within speaker scaling of pitch does not result in a bimodal distribution (as expected). See Figure 2.7.
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_pitch
)
) +
geom_histogram(binwidth=0.1) +
labs(
title = "Within-speaker scaled mean pitch distribution",
x = "Scaled mean pitch"
)
Figure 2.7: Within speaker scaled mean pitch distribution.
We also look at articulation rate.
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_art_rate
)
) +
geom_histogram(binwidth=0.1) +
labs(
title = "Within speaker scaled articulation rate distribution",
x = "Scaled utterance articulation rate"
)
Figure 2.8: Within speaker scaled utterance articulation rate distribution.
Both Figure 2.7, and to a lesser extent, Figure @ref(fig(speaker-scaled-artrate-distribution), show the presence of outliers in the distribution. In both cases, and for consistency, we apply the 2.5 standard deviation rule.
qb_vowels <- qb_vowels %>%
filter(
abs(speaker_scaled_art_rate) <= 2.5,
abs(speaker_scaled_pitch) <= 2.5
)
Some of our models will use a random effect structure to capture the differences between speakers. Consequently, we look at the distributions by speaker.
We know that different speakers contribute radically different amounts to the dataset (Figure 2.9).
qb_vowels %>%
ggplot(
aes(
x = speaker_length
)
) +
geom_histogram(binwidth=100) +
labs(
title = "Distribution of speaker monologue lengths",
x = "Speaker monolgoue length"
)
Figure 2.9: Distribution of speaker monologue lengths.
We can look at the distribution of F1s by speaker to see if any stick out as having not many tokens and, consequently, extreme looking distributions when scaled.
large_plot <- qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_F1,
colour = Speaker
)
) +
facet_wrap(vars(Vowel)) +
geom_freqpoly(stat="density") +
theme(legend.position = "None") +
labs(
title = "Speaker F1 distributions for each vowel",
x = "Scaled F1 value"
)
large_plot
Figure 2.10: Speaker F1 distributions for each vowel.
Figure 2.10 shows some sharp spikes which suggest a lack of data points for certain speakers for certain vowels. We’ll insist on speakers having at least five tokens for each vowel.
qb_vowels <- qb_vowels %>%
group_by(Speaker, Vowel) %>%
mutate(
n_obs = n()
) %>%
ungroup() %>%
group_by(Speaker) %>%
mutate(
min_obs = min(n_obs)
) %>%
filter(
min_obs >= 5
)
Figure 2.11 shows the distribution after filtering. start still has a speaker with a strange distribution, but we have no obvious criterion for removing them.
large_plot <- qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_F1,
colour = Speaker
)
) +
facet_wrap(vars(Vowel)) +
geom_freqpoly(stat="density") +
theme(legend.position = "None") +
labs(
title = "Speaker F1 distributions for each vowel",
subtitle = "After filtering steps",
x = "Scaled F1 value"
)
large_plot
Figure 2.11: Speaker F1 distributions for each vowel (filtered).
In the paper, we compare the extent of amplitude variation in our data with that in laboratory studies. To get a sense how how much variation there is here, we produce some plots. First we look at the variation at the speaker level.
ranges <- qb_vowels %>%
group_by(Speaker) %>%
summarise(
max_amp = max(intensity_max, na.rm=TRUE),
min_amp = min(intensity_max, na.rm=TRUE),
range_amp = max_amp - min_amp
) %>%
ungroup()
summary(ranges$range_amp)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.667 15.800 18.070 18.353 20.431 37.256
qb_vowels %>%
group_by(Speaker) %>%
mutate(
intensity_centered = intensity_max - mean(intensity_max, na.rm=TRUE)
) %>%
ungroup() %>%
ggplot(
aes(
x = intensity_centered,
colour = Speaker,
after_stat(density)
)
) +
geom_freqpoly(binwidth=1) +
labs(
title = "Amplitude variation of Each Speaker (Db)",
) +
theme(
legend.position = "none"
)
qb_vowels %>%
group_by(Speaker) %>%
mutate(
intensity_centered = intensity_max - mean(intensity_max, na.rm=TRUE)
) %>%
ungroup() %>%
ggplot(
aes(
x = intensity_centered,
after_stat(density)
)
) +
geom_freqpoly(binwidth=1) +
labs(
title = "Amplitude Variation in Data (Db)"
) +
theme(
legend.position = "none"
)
And then the amount of variation at the 60 second interval level.
# Now degree of amplitude variation across 60 second intervals.
intervals <- qb_vowels %>%
# Divide up the time variable into 60 second and 240 second intervals.
group_by(Speaker) %>%
mutate(
interval_60 = as.numeric(
as.factor(
cut(
time,
breaks = seq(0, max(time) + 60, 60))
)
)*60
) %>%
# Trim terminal intervals with insufficient data. Replaces bad interval values with NA.
# note: still grouped by Speaker
mutate(
speaker_length = max(time),
remaining_from_start_60 = speaker_length - (interval_60 - 60),
interval_60 = if_else(remaining_from_start_60 >= 45, interval_60, NA_real_),
) %>%
## Create centered intensity
mutate(
intensity_centered = intensity_max - mean(intensity_max, na.rm=TRUE)
) %>%
ungroup() %>%
# Take summary value for formants for 60s intervals.
group_by(Speaker, interval_60) %>%
mutate(
centered_amp_60 = mean(intensity_centered, na.rm=TRUE)
) %>%
ungroup()
intervals %>%
group_by(Speaker, interval_60) %>%
summarise(
centered_amp_60 = first(centered_amp_60)
) %>%
ungroup() %>%
ggplot(
aes(
x = centered_amp_60,
after_stat(density)
)
) +
geom_freqpoly(binwidth=1) +
theme(
legend.position = "none"
)
intervals %>%
group_by(Speaker) %>%
summarise(
max_amp = max(centered_amp_60, na.rm=TRUE),
min_amp = min(centered_amp_60, na.rm=TRUE),
speaker_range = max_amp-min_amp
) %>%
pull(speaker_range) %>%
summary()
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3777 2.8734 4.1153 4.5503 5.6388 21.1940
Before modelling, we want some visual sense of shifts over time in the data.
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_time,
y = speaker_scaled_F1
)
) +
geom_smooth() +
facet_wrap(vars(Vowel)) +
labs(
title = "Changes in F1 over course of monologue",
x = "Time (scaled)",
y = "F1 (scaled)"
)
Figure 2.12: Changes in F1 over course of monologue.
We have some evidence in Figure 2.12 for a decrease in F1 for dress, fleece, foot, lot, nurse, strut, thought, and trap. This may, in turn, be explained by changes in amplitude. Figure 2.13. shows the relationship of amplitude to (scaled) time. We see a small decrease in amplitude over the course of the monologue.
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_time,
y = speaker_scaled_amp_max
)
) +
geom_smooth() +
labs(
title = "Change in amplitude over course of monologue",
x = "Time (scaled)",
y = "Maximum word amplitude (scaled)"
)
Figure 2.13: Change in amplitude over course of monologue.
One final thing to look at here is whether this apparent effect is there for speakers of different lengths. We cut speaker lengths in to 0-10m, 10-20m, and 20+m.
qb_vowels <- qb_vowels %>%
ungroup() %>%
mutate(
speaker_length_fact = cut(
speaker_length,
breaks = c(0, 600, 1200, max(speaker_length)),
labels = c("short (-10m)", "medium (10-20m)", "long (20m+)")
),
speaker_length_fact = fct_relevel(
speaker_length_fact,
c("short (-10m)", "medium (10-20m)", "long (20m+)")
)
)
We look at the number of speakers in each category.
qb_vowels %>%
group_by(speaker_length_fact) %>%
summarise(
n = n_distinct(Speaker)
)
There are no categories with radically low numbers of speakers.
Figure 2.14
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_time,
y = speaker_scaled_amp_max,
colour = speaker_length_fact
)
) +
geom_smooth() +
labs(
title = "Changes in amplitude over time by speaker length",
colour = "Speaker length",
x = "Time (scaled)",
y = "Max word amplitude (scaled)"
)
Figure 2.14: Changes in amplitude over time by speaker length.
All seem to be dropping at the end of the monologue. All seem to cross the overall mean amplitude some time after half way and have a quite precipitous drop at the end.
qb_vowels %>%
ggplot(
aes(
x = speaker_scaled_time,
y = speaker_scaled_pitch
)
) +
geom_smooth() +
labs(
title = "Pitch over the course of monologue.",
y = "Pitch (scaled)",
x = "Time (scaled)"
)
Figure 2.15: Change in pitch over the course of monologue.
We also see a decline in pitch over the course of monologues (Figure 2.15).
We start with some straightforward linear models. First, we use position in the monologue as our predictor.
linear.fit <- lm(speaker_scaled_F1 ~ speaker_scaled_time, data=qb_vowels)
summary(linear.fit)
##
## Call:
## lm(formula = speaker_scaled_F1 ~ speaker_scaled_time, data = qb_vowels)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2825 -0.6506 -0.0079 0.6413 4.2431
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.040975 0.005883 6.965 3.3e-12 ***
## speaker_scaled_time -0.104655 0.010307 -10.153 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9784 on 109341 degrees of freedom
## Multiple R-squared: 0.000942, Adjusted R-squared: 0.0009328
## F-statistic: 103.1 on 1 and 109341 DF, p-value: < 2.2e-16
A simple model with scaled time as predictor strongly indicates a reduction in F1 over time. However, with an adjusted R-squared of 0.00059, this model explains almost none of the variation in F1 in the dataset.
We will get a little bit more sophisticated before engaging in model diagnostics.
linear.fit.2 <- lm(
speaker_scaled_F1 ~
speaker_scaled_time*scaled_length +
speaker_scaled_art_rate +
speaker_scaled_amp_max +
speaker_scaled_pitch,
data=qb_vowels
)
summary(linear.fit.2)
##
## Call:
## lm(formula = speaker_scaled_F1 ~ speaker_scaled_time * scaled_length +
## speaker_scaled_art_rate + speaker_scaled_amp_max + speaker_scaled_pitch,
## data = qb_vowels)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1982 -0.6430 -0.0197 0.6239 4.3273
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.011508 0.005869 1.961 0.049912 *
## speaker_scaled_time -0.056796 0.010271 -5.530 3.22e-08 ***
## scaled_length 0.018010 0.005738 3.139 0.001698 **
## speaker_scaled_art_rate -0.013432 0.003089 -4.349 1.37e-05 ***
## speaker_scaled_amp_max 0.176110 0.003756 46.889 < 2e-16 ***
## speaker_scaled_pitch -0.014982 0.004250 -3.525 0.000424 ***
## speaker_scaled_time:scaled_length -0.037333 0.010082 -3.703 0.000213 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9663 on 109336 degrees of freedom
## Multiple R-squared: 0.02553, Adjusted R-squared: 0.02547
## F-statistic: 477.3 on 6 and 109336 DF, p-value: < 2.2e-16
Our second model includes the across-speaker scaled length of the speaker and interacts it with time scaled. The idea here is that longer speakers may exhibit a more extreme effect of time on their F1, particularly if fatigue is part of the story. We also include articulation rate and amplitude taken at midpoint.
The largest factor here is, by a multiple of 10, amplitude. Scaled time, by itself, is not significant, but the interaction with speaker length does have an effect. This is not surprising, insofar as this is a continuous by continuous interaction where the baseline is a speaker length of 0. The effect of scaled time in reducing F1 seems to increase as speakers increase in length. The articulation rate term suggests that as articulation rate increases, the F1 will decrease, and similar for the pitch.
This model moves us from an R squared of 0.00059 to 0.019. This is a major increase in descriptive power.
We check the diagnostic plots.
plot(linear.fit.2)
There is no extreme departure from normality of residuals and no data points
with excessive leverage. In the tails, the error distribution is not entirely
normal, but is within acceptable limits.
It is easier to interpret a speaker length factor than the continuous by continuous interaction, so we refit the model.
linear.fit.2.fact <- lm(
scaled_F1 ~
speaker_scaled_time*speaker_length_fact +
speaker_scaled_art_rate +
speaker_scaled_amp_max +
speaker_scaled_pitch,
data=qb_vowels
)
summary(linear.fit.2.fact)
##
## Call:
## lm(formula = scaled_F1 ~ speaker_scaled_time * speaker_length_fact +
## speaker_scaled_art_rate + speaker_scaled_amp_max + speaker_scaled_pitch,
## data = qb_vowels)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.76046 -0.71915 -0.01502 0.76985 2.96135
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -0.011465 0.011571
## speaker_scaled_time 0.042369 0.020270
## speaker_length_factmedium (10-20m) -0.005132 0.015140
## speaker_length_factlong (20m+) 0.013291 0.015619
## speaker_scaled_art_rate -0.007767 0.003237
## speaker_scaled_amp_max -0.058774 0.003936
## speaker_scaled_pitch 0.060664 0.004454
## speaker_scaled_time:speaker_length_factmedium (10-20m) -0.013882 0.026535
## speaker_scaled_time:speaker_length_factlong (20m+) 0.057375 0.027357
## t value Pr(>|t|)
## (Intercept) -0.991 0.3218
## speaker_scaled_time 2.090 0.0366 *
## speaker_length_factmedium (10-20m) -0.339 0.7346
## speaker_length_factlong (20m+) 0.851 0.3948
## speaker_scaled_art_rate -2.400 0.0164 *
## speaker_scaled_amp_max -14.933 <2e-16 ***
## speaker_scaled_pitch 13.621 <2e-16 ***
## speaker_scaled_time:speaker_length_factmedium (10-20m) -0.523 0.6009
## speaker_scaled_time:speaker_length_factlong (20m+) 2.097 0.0360 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.013 on 109334 degrees of freedom
## Multiple R-squared: 0.003498, Adjusted R-squared: 0.003425
## F-statistic: 47.97 on 8 and 109334 DF, p-value: < 2.2e-16
We see that longer speakers have a stronger time scaled effect than shorter speakers. We also see that amplitude remains a major factor. There is no evidence here for an effect of F1 drop over time distinct from the drop in amplitude for all but the longest speakers.
Before incorporating by-speaker random effects, it is worth looking at models
for each vowel to determine whether the same predictors come out as significant.
We do this by using the purrr method of nesting, fitting a model to each
nested dataframe, and extracting model summary information.
vowel_linear_models <- qb_vowels %>%
# Group by vowel and nest to create a column of dataframes corresponding
# to each vowel.
group_by(Vowel) %>%
nest() %>%
# Apply the linear model (same structure as linear.fit.2)
mutate(
model = map(
data,
~ lm(speaker_scaled_F1 ~
speaker_scaled_time*scaled_length +
speaker_scaled_art_rate +
speaker_scaled_amp_max +
speaker_scaled_pitch,
data = .x)
),
# Extract the coefficients for each variable
coefficients = map(model, tidy),
# List the variables with p-values less than or equal to 0.05.
significant_variables = map(coefficients, ~ .x %>% filter(p.value <= 0.05))
) %>%
# Select the significant variables and 'unnest' so that each vowel has a row
# for each significant variable.
select(
Vowel, significant_variables
) %>%
unnest(significant_variables)
We then output the models where speaker_scaled_time or the interaction
speaker_scaled_time:scaled_length come out as statistically significant at
the 0.05 level.
# Output the models where scaled_time or the interaction with speaker
# length is significant.
vowel_linear_models %>%
filter(
term %in% c("speaker_scaled_time", "speaker_scaled_time:scaled_length")
)
The models for kit, trap, and nurse, thought, and lot suggest that something is happening with scaled_time distinct from changes in amplitude. This is particularly interesting given that both kit and start show almost no movement in average F1 over time (Figure 2.12).
We now output all the models where scaled_amp_max appears as a significant
predictor.
vowel_linear_models %>%
filter(term == "speaker_scaled_amp_max")
Amplitude is taken to be a significant predictor by all vowels. Note also, that the effect is in the same direction for each, with lot and strut having the strongest effects by magnitude.
vowel_linear_models %>%
filter(term == "speaker_scaled_pitch")
The pitch seems to also be associated with F1 for dress, nurse, kit, strut, lot, and foot, with the strongest effect for strut.
Finally, we look at articulation rate:
vowel_linear_models %>%
filter(term == "speaker_scaled_art_rate")
This predictor is taken to be significant for eight monophthongs. But note that the direction of the effect is different for goose and thought. This is consistent with a decreased vowel space when articulation rate is high (with goose and thought lowering while the other vowels rise.
We now move to linear models which incorporate random effects structures. This allows us to capture the fact that our data comes from different speakers who, we may assume, are drawn from a roughly normal distribution of possible speakers.
The linear models suggests that there is variation in the effects of our predictors by each vowel. Vowel levels are repeatable, so we do not include them as random effects. We only include speaker as a random effect.
We model using raw F1 frequency, as this makes the differences between the vowels easier for the model to detect. We will use random intercepts for each speaker-vowel combination.
NB: to refit this model yourself, change eval to TRUE in the following code
chunk options.
lmer_fit <- lmer(
F1_50 ~
Vowel +
speaker_scaled_time * scaled_length +
speaker_scaled_amp_max * Vowel +
speaker_scaled_art_rate * Vowel +
speaker_scaled_pitch * Vowel +
participant_gender +
(1 + Vowel|Speaker),
data = qb_vowels)
# save model
write_rds(lmer_fit, here('models', 'lmer_fit.rds'))
We load a model we have already fit (or, if the previous cell has been run, the model which has just been saved).
# load model
lmer_fit <- read_rds(here('models', 'lmer_fit.rds'))
The easiest way to look at the model effects is to look at the output of the
anova function. This outputs test values for variable terms rather than for
each individual level of our Vowel factor.
anova(lmer_fit)
The majority of the variation handled by the model’s fixed effects comes from
the different vowels. The F value column suggests that amplitude, scaled by
speaker, is significant. It, and the interaction with Vowel explain a very
large proportion of the variation in the data which is explained by the model.
We can also look at the detailed output.
summary(lmer_fit)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## F1_50 ~ Vowel + speaker_scaled_time * scaled_length + speaker_scaled_amp_max *
## Vowel + speaker_scaled_art_rate * Vowel + speaker_scaled_pitch *
## Vowel + participant_gender + (1 + Vowel | Speaker)
## Data: qb_vowels
##
## REML criterion at convergence: 852275.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.2332 -0.5744 -0.0179 0.5743 6.2984
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Speaker (Intercept) 1456.7 38.17
## VowelFLEECE 1680.2 40.99 -0.60
## VowelFOOT 1406.4 37.50 -0.60 0.78
## VowelGOOSE 1395.2 37.35 -0.62 0.90 0.82
## VowelKIT 2218.6 47.10 -0.49 0.78 0.87 0.73
## VowelLOT 1898.3 43.57 -0.38 0.46 0.65 0.44 0.63
## VowelNURSE 627.8 25.05 -0.08 0.23 0.09 0.13 0.00 -0.04
## VowelSTART 7803.6 88.34 -0.36 0.51 0.54 0.46 0.68 0.63
## VowelSTRUT 5364.0 73.24 -0.29 0.47 0.60 0.41 0.72 0.70
## VowelTHOUGHT 601.5 24.53 -0.48 0.51 0.66 0.57 0.43 0.46
## VowelTRAP 1189.1 34.48 -0.02 0.08 -0.01 -0.07 0.22 0.38
## Residual 3140.1 56.04
##
##
##
##
##
##
##
##
## -0.12
## -0.11 0.88
## 0.50 0.25 0.28
## 0.33 0.45 0.44 0.10
##
## Number of obs: 77901, groups: Speaker, 171
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 421.99564 3.44054 122.654
## VowelFLEECE -8.61423 3.24913 -2.651
## VowelFOOT 65.32696 3.25347 20.079
## VowelGOOSE -21.27759 3.06497 -6.942
## VowelKIT 87.46343 3.70494 23.607
## VowelLOT 168.25271 3.47590 48.406
## VowelNURSE 12.12402 2.29488 5.283
## VowelSTART 366.89692 6.88447 53.293
## VowelSTRUT 265.08340 5.67146 46.740
## VowelTHOUGHT 10.34324 2.12530 4.867
## VowelTRAP 149.05113 2.80505 53.137
## speaker_scaled_time -1.45227 0.73004 -1.989
## scaled_length -1.24112 2.90666 -0.427
## speaker_scaled_amp_max 2.46245 0.64666 3.808
## speaker_scaled_art_rate 0.09395 0.50336 0.187
## speaker_scaled_pitch 1.80682 0.73226 2.467
## participant_genderM -33.21261 4.81888 -6.892
## speaker_scaled_time:scaled_length -1.29869 0.70726 -1.836
## VowelFLEECE:speaker_scaled_amp_max 1.04060 0.98074 1.061
## VowelFOOT:speaker_scaled_amp_max 7.54225 1.74881 4.313
## VowelGOOSE:speaker_scaled_amp_max 5.36673 1.25829 4.265
## VowelKIT:speaker_scaled_amp_max 5.54461 0.95933 5.780
## VowelLOT:speaker_scaled_amp_max 19.24589 1.05992 18.158
## VowelNURSE:speaker_scaled_amp_max 1.94509 1.48429 1.310
## VowelSTART:speaker_scaled_amp_max 7.41648 1.32731 5.588
## VowelSTRUT:speaker_scaled_amp_max 19.47301 0.98573 19.755
## VowelTHOUGHT:speaker_scaled_amp_max 1.77404 1.18602 1.496
## VowelTRAP:speaker_scaled_amp_max 4.63630 1.08204 4.285
## VowelFLEECE:speaker_scaled_art_rate 0.05099 0.77399 0.066
## VowelFOOT:speaker_scaled_art_rate -3.08327 1.44553 -2.133
## VowelGOOSE:speaker_scaled_art_rate 2.68636 1.02888 2.611
## VowelKIT:speaker_scaled_art_rate -2.79113 0.77505 -3.601
## VowelLOT:speaker_scaled_art_rate -3.17524 0.87903 -3.612
## VowelNURSE:speaker_scaled_art_rate 0.21605 1.14143 0.189
## VowelSTART:speaker_scaled_art_rate -3.48098 1.05748 -3.292
## VowelSTRUT:speaker_scaled_art_rate -2.76832 0.80442 -3.441
## VowelTHOUGHT:speaker_scaled_art_rate 1.20752 0.93501 1.291
## VowelTRAP:speaker_scaled_art_rate -3.38389 0.86798 -3.899
## VowelFLEECE:speaker_scaled_pitch -3.24552 1.11643 -2.907
## VowelFOOT:speaker_scaled_pitch -7.10082 1.91541 -3.707
## VowelGOOSE:speaker_scaled_pitch -1.61644 1.48691 -1.087
## VowelKIT:speaker_scaled_pitch -5.70438 1.06770 -5.343
## VowelLOT:speaker_scaled_pitch -8.28547 1.22217 -6.779
## VowelNURSE:speaker_scaled_pitch 0.72154 1.59803 0.452
## VowelSTART:speaker_scaled_pitch -1.07707 1.54581 -0.697
## VowelSTRUT:speaker_scaled_pitch -12.23869 1.16369 -10.517
## VowelTHOUGHT:speaker_scaled_pitch -1.07150 1.33284 -0.804
## VowelTRAP:speaker_scaled_pitch -0.01169 1.24410 -0.009
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00420057 (tol = 0.002, component 1)
Note that model’s convergence is somewhat higher (0.0042) than the default tolerance (0.002), but is reasonably close. As this isn’t the main model we will use in the paper, we will not engage in any further attempts to get convergence within the threshold.
Breaking up the results by vowel reveals some important patterns. In particular,
we see that lot and strut dominate when we look at the
scaled_amp_max effect. The effect of scaled_amp_max on F1 is positive apart
from the very low magnitude negative coefficient for START.
We check the diagnostic plots, as above.
qqnorm(resid(lmer_fit))
qqline(resid(lmer_fit))
The model is struggling at the tails. This may be corrected by allowing for
non-linear relationships when we turn to GAMMs.
plot(lmer_fit)
We know that amplitude, pitch, and articulation rate have relationships to one another. But are these causing problems for our estimated coefficients. We look at the variance inflation factor.
vif(lmer_fit)
## GVIF Df GVIF^(1/(2*Df))
## Vowel 1.184820 10 1.008516
## speaker_scaled_time 1.054574 1 1.026925
## scaled_length 1.034398 1 1.017053
## speaker_scaled_amp_max 8.035549 1 2.834704
## speaker_scaled_art_rate 5.565775 1 2.359189
## speaker_scaled_pitch 7.848042 1 2.801436
## participant_gender 1.015946 1 1.007942
## speaker_scaled_time:scaled_length 1.061346 1 1.030217
## Vowel:speaker_scaled_amp_max 108.863636 10 1.264283
## Vowel:speaker_scaled_art_rate 6.238248 10 1.095855
## Vowel:speaker_scaled_pitch 106.417448 10 1.262847
A GVIF above 2.24 represents some multicollinearity, a GVIF above 3.16 indicates high colinearity. The variance inflation here may come from the correlation between amplitude and pitch (0.2758935). Since the GVIF is not extreme, we will leave all variable in.
We can now plot the amplitude and F1 relationship for each vowel according for the generalised linear model.
# Define new data to generate model predictions
new_data <- tibble(
speaker_scaled_pitch = 0,
speaker_scaled_art_rate = 0,
scaled_length = 0,
Vowel = rep(vowels, each = 1000),
participant_gender = "F",
# Assuming we are in the middle of the monologue
speaker_scaled_time = 0.5,
speaker_scaled_amp_max = rep(seq(-3, 3, length.out = 1000), times=11)
)
new_data <- new_data %>%
mutate(
prediction = predict(lmer_fit, newdata=new_data, re.form=NA),
)
new_data %>%
ggplot(
aes(
x = speaker_scaled_amp_max,
y = prediction,
colour = Vowel
)
) +
geom_line() +
labs(
title = "F1 and Amplitude by Vowel",
subtitle = glue(
"Prediction, controlling for speaker gender, length, time, articulation rate,",
"and pitch"
),
colour = "Speaker length",
y = "Predicted F1 (hz)",
x = "Maximum amplitude (scaled)"
) +
scale_colour_manual(
values = vowel_colours_with_foot
)
Figure 2.16: Relationship between F1 and amplitude by vowel.
Figure 2.16 shows predicted F1 values for different vowels at different maximum amplitudes. We see that the effect, measured in terms of raw frequency is very large for strut and lot. What is also interesting is that we see some overlap between vowels, so that, for instance, at low amplitude LOT has a lower F1 than trap, but at moderate to high volumes trap has a lower F1 than LOT. It also appears that, at as the amplitude increases, the F1 values of goose and fleece will get closer and closer together. At low amplitude, the F1 value of thought, nurse, fleece and dress coincide, but become distinguished from each other at higher amplitudes. On the other hand, at low amplitudes we find GOOSE more distinguished in height from the other high vowels, but coinciding with fleece at high amplitudes.
We now turn to GAMMs, a more flexible approach to modelling which is particularly effective at capturing non-linear relationships. The failure to handle tails of our error distribution with linear mixed models suggests that we might benefit from allowing non-linearities in to our models.
First, we fix the representation of the data for the mgcv package, which
requires factors to be represented as such rather than interpreting character
vectors as factors.
qb_vowels <- qb_vowels %>%
mutate(
Vowel = as.factor(Vowel),
Speaker = as.factor(Speaker),
participant_gender = as.factor(participant_gender)
)
We begin by exploring structures using fREML to fit our models. We then use ML in order to carry out model-comparison based significance testing. The full model we fit using ML is the model from which we generate the effect plots in the paper.
We are not performing significance tests for distinctions between vowels, so we do not fit difference smooths (as in Sókuthy).
We fit a series of models here using fREML. These can be fit quickly, however, they are not suitable for model comparison-based testing when fixed effects are varied (Sóskuthy 2021, 8). In order to do this, we will fit a series of models with ML. These will take much longer to fit.
Our first structure is used in the main paper for vowel space visualisations in which the effect of amplitude and articulation rate are compared.
Our second structure, which adds random smooths over the time variable for each speaker, is given as part of the process we carried out to explore potential model structures. It is not used in the main paper. We did not find that it significantly improved performance for the added model complexity.
Our third structure explores the possibility of fitting the model using a scaled t-distribution as our error model rather than a Gaussian. This is found to improve the behaviour of our residuals, but it does not significantly affect the interpretation of our results. Moreover, attempts to fit the model with scaled t residuals were unsuccessful with ML due to limitations in computational power.
We begin by following the same structure as the previous mixed model.
NB: to refit this model yourself, change eval to TRUE in the following code
chunk options.
gamm_fit <- bam(
F1_50 ~
participant_gender +
Vowel +
s(speaker_scaled_time, by=Vowel) +
s(speaker_scaled_art_rate, by=Vowel) +
s(speaker_scaled_amp_max, by=Vowel) +
s(speaker_scaled_pitch, by=Vowel) +
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="fREML",
discrete = TRUE,
nthreads = 8
)
write_rds(gamm_fit, here('models', 'gamm_fit_s1.rds'))
# Also takes a long time to calculate summary.
# We calculate it and save it now.
gamm_fit_summary <- summary(gamm_fit)
write_rds(gamm_fit_summary, here('models', 'gamm_fit_s1_summary.rds'))
We load the pre-calculated/fit model and summary.
gamm_fit <- read_rds(here('models', 'gamm_fit_s1.rds'))
gamm_fit_summary <- read_rds(here('models', 'gamm_fit_s1_summary.rds'))
gamm_fit_summary
##
## Family: gaussian
## Link function: identity
##
## Formula:
## F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(Speaker,
## bs = "re") + s(Speaker, Vowel, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 434.637 4.115 105.617 < 2e-16 ***
## participant_genderM -55.282 5.637 -9.807 < 2e-16 ***
## VowelFLEECE -13.294 3.983 -3.338 0.000844 ***
## VowelFOOT 62.689 4.477 14.001 < 2e-16 ***
## VowelGOOSE -21.115 4.193 -5.036 4.76e-07 ***
## VowelKIT 84.363 3.981 21.190 < 2e-16 ***
## VowelLOT 155.426 4.098 37.928 < 2e-16 ***
## VowelNURSE 13.788 4.197 3.285 0.001018 **
## VowelSTART 355.660 4.002 88.862 < 2e-16 ***
## VowelSTRUT 251.938 4.102 61.414 < 2e-16 ***
## VowelTHOUGHT 10.656 4.066 2.621 0.008775 **
## VowelTRAP 144.036 4.136 34.829 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(speaker_scaled_time):VowelDRESS 1.807 2.257 3.772 0.018803 *
## s(speaker_scaled_time):VowelFLEECE 1.003 1.005 2.740 0.097453 .
## s(speaker_scaled_time):VowelFOOT 1.003 1.005 0.996 0.319046
## s(speaker_scaled_time):VowelGOOSE 1.001 1.002 0.000 0.996387
## s(speaker_scaled_time):VowelKIT 5.733 6.889 2.647 0.012835 *
## s(speaker_scaled_time):VowelLOT 2.549 3.177 9.334 3.17e-06 ***
## s(speaker_scaled_time):VowelNURSE 1.001 1.002 0.537 0.463580
## s(speaker_scaled_time):VowelSTART 1.003 1.005 0.017 0.904810
## s(speaker_scaled_time):VowelSTRUT 6.321 7.472 4.470 0.000111 ***
## s(speaker_scaled_time):VowelTHOUGHT 1.002 1.004 6.115 0.013364 *
## s(speaker_scaled_time):VowelTRAP 2.489 3.104 6.832 0.000110 ***
## s(speaker_scaled_art_rate):VowelDRESS 1.002 1.003 0.032 0.861188
## s(speaker_scaled_art_rate):VowelFLEECE 1.164 1.310 1.506 0.171256
## s(speaker_scaled_art_rate):VowelFOOT 1.001 1.003 7.013 0.008075 **
## s(speaker_scaled_art_rate):VowelGOOSE 1.001 1.003 24.048 8.40e-07 ***
## s(speaker_scaled_art_rate):VowelKIT 1.001 1.003 40.826 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelLOT 1.002 1.005 25.488 4.62e-07 ***
## s(speaker_scaled_art_rate):VowelNURSE 1.001 1.003 0.161 0.689283
## s(speaker_scaled_art_rate):VowelSTART 1.001 1.002 18.858 1.39e-05 ***
## s(speaker_scaled_art_rate):VowelSTRUT 3.587 4.475 6.527 1.42e-05 ***
## s(speaker_scaled_art_rate):VowelTHOUGHT 1.001 1.002 5.908 0.015051 *
## s(speaker_scaled_art_rate):VowelTRAP 1.002 1.003 11.084 0.000860 ***
## s(speaker_scaled_amp_max):VowelDRESS 1.001 1.002 24.459 7.31e-07 ***
## s(speaker_scaled_amp_max):VowelFLEECE 2.351 2.995 20.272 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelFOOT 1.041 1.080 93.897 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelGOOSE 1.775 2.248 43.811 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelKIT 3.685 4.605 55.138 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelLOT 3.889 4.843 270.462 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelNURSE 1.001 1.002 22.382 2.23e-06 ***
## s(speaker_scaled_amp_max):VowelSTART 1.224 1.416 71.417 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelSTRUT 5.481 6.600 223.904 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelTHOUGHT 1.404 1.714 19.053 1.81e-06 ***
## s(speaker_scaled_amp_max):VowelTRAP 3.669 4.589 45.878 < 2e-16 ***
## s(speaker_scaled_pitch):VowelDRESS 3.863 4.804 5.465 7.72e-05 ***
## s(speaker_scaled_pitch):VowelFLEECE 2.520 3.197 3.902 0.008265 **
## s(speaker_scaled_pitch):VowelFOOT 3.246 4.056 8.454 9.65e-08 ***
## s(speaker_scaled_pitch):VowelGOOSE 3.033 3.817 3.604 0.006771 **
## s(speaker_scaled_pitch):VowelKIT 5.792 6.910 25.970 < 2e-16 ***
## s(speaker_scaled_pitch):VowelLOT 4.788 5.855 23.215 < 2e-16 ***
## s(speaker_scaled_pitch):VowelNURSE 3.175 3.979 5.173 0.000404 ***
## s(speaker_scaled_pitch):VowelSTART 1.002 1.003 3.106 0.077653 .
## s(speaker_scaled_pitch):VowelSTRUT 7.884 8.638 19.806 < 2e-16 ***
## s(speaker_scaled_pitch):VowelTHOUGHT 3.072 3.864 5.081 0.000563 ***
## s(speaker_scaled_pitch):VowelTRAP 3.795 4.725 3.587 0.003102 **
## s(Speaker) 192.183 214.000 2230.801 < 2e-16 ***
## s(Vowel,Speaker) 1913.920 2361.000 37.150 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.759 Deviance explained = 76.4%
## fREML = 6.2453e+05 Scale est. = 5061.3 n = 109343
The summary output here has a great deal of detail. We will focus in on specific
aspects of the model as we go. At this point, it is worth noting that all of
the terms here are significant, at least for one level of our Vowel factor.
The model here first a smooth for each of our predictors and then enables a
difference from this smooth for each vowel.
We will look at amplitude first, as it is our main variable of interest. The code below can be modified to get and plot predictions for the other variables.
gamm_fit_preds <- get_predictions(
gamm_fit,
cond = list(
'speaker_scaled_amp_max' = seq(-2.5, 2.5, length.out=200),
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
gamm_fit_preds %>%
mutate(# Enable plot to distinguish between high and low vowels (for faceting)
height = if_else(Vowel %in% high_vowels, "high", "low or mid"),
Vowel = factor(Vowel, levels = vowels)
) %>%
ggplot(
aes(
y = fit,
x = speaker_scaled_amp_max,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.3,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = "F1 and Amplitude by Vowel",
x = "Scaled amplitude",
y = "F1"
)
Figure 2.17: F1 and amplitude by vowel.
# facet_wrap(
# facets=vars(height), scales = "free"
# )
Figure 2.17 looks very similar to Figure 2.16. We have the crossing of lot and trap, for instance, and, broadly speaking, increased F1 with increased relative amplitude.
We now run diagnostic checks on the model:Figure 2.18: GAM model checks.
A few things are worth noting here. First, we don’t have any indication that
we have our \(k\) values set too low. The check for this can be found at the bottom
of the text output of gam.check above, where we need both an edf value close
to k' and a p-value less than 1. This is not the case for any of our variables.
NB: \(k\) can be thought of as an upper limit on ‘wiggliness’. If k is too low,
then we are not giving the GAM enough freedom to fit the true relationship.
The QQ plot again suggests that we are struggling at the tails of our distribution of residuals. This is not particularly worrying, although we will investigate an alternative model for our residuals below (scaled-t rather than Gaussian). The response vs. fitted variables plot looks fine. As expected, given the QQ plot, the histogram of residuals has quite heavy tails.
An important factor in gamm modelling is to ensure that the model specification
enables sufficient ‘wiggliness’. An indication that this is not the case is that
estimated degrees of freedom values are close to the ‘k’ value. All of our
smooths take the default number of knows (k=10). We look at the text output
of gam.check to ensure that the estimated degrees of freedom are lower than
our k value. This is indeed the case.
Note the absence of autocorrelation in the model residuals (Figure 2.19).
acf_resid(gamm_fit, split_pred = 'Speaker')
Figure 2.19: Autocorrelation in model residuals.
We will take this to apply to the other models of the data we fit.
We can also set up a model with by-speaker random smooths. This takes a very long time to run. By-speaker smooths are preferable in so far as they enable the model to more accurately capture how much independent information we have. Our previous model structure does not understand our datapoints as occurring within an ongoing monologue for a particular speaker.
NB: to refit this model yourself, change eval to TRUE in the following code
chunk options.
gamm_fit_by_speaker <- bam(
F1_50 ~
participant_gender +
Vowel +
s(speaker_scaled_time, by = Vowel) +
s(speaker_scaled_art_rate, by=Vowel) +
s(speaker_scaled_amp_max, by=Vowel) +
s(speaker_scaled_pitch, by=Vowel) +
s(speaker_scaled_time, Speaker, bs="fs", k=5, m=1) + #k=5 to reduce load.
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="fREML",
discrete = TRUE,
nthreads = 8
)
write_rds(gamm_fit_by_speaker, here('models', 'gamm_fit_by_speaker.rds'))
# We again calculate and save the summary.
gamm_fit_by_speaker_summary <- summary(gamm_fit_by_speaker)
write_rds(
gamm_fit_by_speaker_summary,
here('models', 'gamm_fit_by_speaker_summary.rds')
)
This cell loads a previously generated version of the model.
gamm_fit_by_speaker <- read_rds(here('models', 'gamm_fit_by_speaker.rds'))
gamm_fit_by_speaker_summary <- read_rds(
here('models', 'gamm_fit_by_speaker_summary.rds')
)
gamm_fit_by_speaker_summary
##
## Family: gaussian
## Link function: identity
##
## Formula:
## F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(speaker_scaled_time,
## Speaker, bs = "fs", k = 5, m = 1) + s(Speaker, bs = "re") +
## s(Speaker, Vowel, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 434.590 4.112 105.701 < 2e-16 ***
## participant_genderM -55.161 5.632 -9.794 < 2e-16 ***
## VowelFLEECE -13.202 3.970 -3.325 0.000883 ***
## VowelFOOT 62.305 4.471 13.937 < 2e-16 ***
## VowelGOOSE -21.252 4.193 -5.068 4.03e-07 ***
## VowelKIT 84.194 3.978 21.166 < 2e-16 ***
## VowelLOT 155.533 4.094 37.993 < 2e-16 ***
## VowelNURSE 13.481 4.191 3.217 0.001297 **
## VowelSTART 355.197 3.972 89.433 < 2e-16 ***
## VowelSTRUT 252.047 4.098 61.511 < 2e-16 ***
## VowelTHOUGHT 10.440 4.068 2.566 0.010289 *
## VowelTRAP 144.209 4.130 34.916 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(speaker_scaled_time):VowelDRESS 1.616 1.987 2.527 0.068312 .
## s(speaker_scaled_time):VowelFLEECE 1.006 1.011 1.754 0.186140
## s(speaker_scaled_time):VowelFOOT 1.002 1.004 1.116 0.290215
## s(speaker_scaled_time):VowelGOOSE 1.001 1.002 0.005 0.948852
## s(speaker_scaled_time):VowelKIT 5.863 7.020 2.585 0.011066 *
## s(speaker_scaled_time):VowelLOT 2.250 2.798 8.608 4.24e-05 ***
## s(speaker_scaled_time):VowelNURSE 1.001 1.003 0.731 0.392396
## s(speaker_scaled_time):VowelSTART 1.005 1.010 0.019 0.903865
## s(speaker_scaled_time):VowelSTRUT 6.256 7.408 3.886 0.000154 ***
## s(speaker_scaled_time):VowelTHOUGHT 1.001 1.002 5.608 0.017846 *
## s(speaker_scaled_time):VowelTRAP 2.255 2.806 6.238 0.000826 ***
## s(speaker_scaled_art_rate):VowelDRESS 1.001 1.002 0.084 0.773479
## s(speaker_scaled_art_rate):VowelFLEECE 1.005 1.011 2.527 0.110218
## s(speaker_scaled_art_rate):VowelFOOT 1.001 1.001 6.066 0.013765 *
## s(speaker_scaled_art_rate):VowelGOOSE 1.003 1.007 23.176 1.44e-06 ***
## s(speaker_scaled_art_rate):VowelKIT 1.002 1.004 40.579 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelLOT 1.003 1.006 26.385 2.08e-07 ***
## s(speaker_scaled_art_rate):VowelNURSE 1.001 1.002 0.163 0.686753
## s(speaker_scaled_art_rate):VowelSTART 1.002 1.003 18.240 1.91e-05 ***
## s(speaker_scaled_art_rate):VowelSTRUT 3.520 4.393 6.918 8.48e-06 ***
## s(speaker_scaled_art_rate):VowelTHOUGHT 1.001 1.003 6.039 0.013972 *
## s(speaker_scaled_art_rate):VowelTRAP 1.003 1.006 10.909 0.000932 ***
## s(speaker_scaled_amp_max):VowelDRESS 1.001 1.002 24.840 6.43e-07 ***
## s(speaker_scaled_amp_max):VowelFLEECE 2.279 2.904 21.227 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelFOOT 1.019 1.037 96.964 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelGOOSE 1.913 2.431 40.792 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelKIT 3.733 4.661 55.259 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelLOT 3.886 4.838 270.720 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelNURSE 1.002 1.004 21.749 3.20e-06 ***
## s(speaker_scaled_amp_max):VowelSTART 1.005 1.011 91.804 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelSTRUT 5.476 6.593 225.431 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelTHOUGHT 1.558 1.947 16.333 1.09e-06 ***
## s(speaker_scaled_amp_max):VowelTRAP 3.603 4.510 46.685 < 2e-16 ***
## s(speaker_scaled_pitch):VowelDRESS 3.838 4.773 5.322 0.000110 ***
## s(speaker_scaled_pitch):VowelFLEECE 2.398 3.044 4.887 0.002143 **
## s(speaker_scaled_pitch):VowelFOOT 3.277 4.091 8.672 3.81e-07 ***
## s(speaker_scaled_pitch):VowelGOOSE 3.033 3.816 3.601 0.006810 **
## s(speaker_scaled_pitch):VowelKIT 5.813 6.929 26.226 < 2e-16 ***
## s(speaker_scaled_pitch):VowelLOT 4.841 5.910 23.817 < 2e-16 ***
## s(speaker_scaled_pitch):VowelNURSE 3.120 3.911 4.912 0.000753 ***
## s(speaker_scaled_pitch):VowelSTART 1.003 1.007 2.436 0.117712
## s(speaker_scaled_pitch):VowelSTRUT 7.860 8.623 20.369 < 2e-16 ***
## s(speaker_scaled_pitch):VowelTHOUGHT 3.064 3.853 4.974 0.000679 ***
## s(speaker_scaled_pitch):VowelTRAP 3.838 4.773 3.666 0.002522 **
## s(speaker_scaled_time,Speaker) 425.051 1078.000 452.834 0.022190 *
## s(Speaker) 96.083 214.000 0.897 < 2e-16 ***
## s(Vowel,Speaker) 1914.372 2361.000 10.702 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.762 Deviance explained = 76.7%
## fREML = 6.2429e+05 Scale est. = 5015.6 n = 109343
We again visualise:
gamm_fit_preds_by_speaker <- get_predictions(
gamm_fit_by_speaker,
cond = list(
'speaker_scaled_amp_max' = seq(-2.5, 2.5, length.out=200),
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.488817037310488.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): -0.0237747342657432.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): -0.0723808600211757.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker),s(Vowel,Speaker)
##
gamm_fit_preds_by_speaker %>%
mutate(# Enable plot to distinguish between high and low vowels (for faceting)
height = if_else(Vowel %in% high_vowels, "high", "low or mid")
) %>%
ggplot(
aes(
y = fit,
x = speaker_scaled_amp_max,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.3,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = "F1 and Amplitude by Vowel (by Speaker Smooths Over Time)",
x = "Scaled amplitude",
y = "F1"
)
Figure 2.20: F1 and amplitude by vowel with by by-speaker smooths over time.
# facet_wrap(
# facets=vars(height), scales = "free"
# )
We again run model diagnostics:
gam.check(gamm_fit_by_speaker)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] -7.611567e-11 -7.831039e-04 -6.873385e-04 -4.269471e-04 -3.232969e-13
## [6] -3.390954e-12 -5.671226e-04 -8.105723e-04 -1.084022e-12 -3.951028e-04
## [11] -6.582179e-12 -4.268056e-04 -3.181207e-04 -3.193595e-04 -7.605391e-04
## [16] -7.023723e-04 -7.369885e-04 -3.502671e-04 -7.385356e-04 -3.812506e-12
## [21] -4.372413e-04 -4.382093e-04 -4.872777e-04 3.103073e-12 -4.900668e-04
## [26] -4.309331e-13 -2.943423e-12 8.684164e-13 -7.964819e-04 -8.866592e-06
## [31] 5.440093e-13 -5.619505e-12 5.160317e-13 -2.529088e-12 -2.448264e-12
## [36] 1.456960e-08 -1.849632e-13 -8.536727e-12 -3.735900e-12 -1.574296e-12
## [41] -4.707181e-04 -1.376144e-11 -2.481571e-12 -7.627232e-13 1.290346e-11
## [46] 2.943068e-11 1.373479e-11 3.519176e-09 -3.859896e-08
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## 1.627120e-01 1.295647e-04 -1.870680e-05 1.129553e-05 1.857063e-03
## 1.295647e-04 7.847806e-04 1.444489e-07 -1.054191e-07 -1.595943e-05
## -1.870680e-05 1.444489e-07 6.868904e-04 1.256478e-08 1.660512e-06
## 1.129553e-05 -1.054191e-07 1.256478e-08 4.267081e-04 -1.182149e-06
## 1.857063e-03 -1.595943e-05 1.660512e-06 -1.182149e-06 7.541779e-01
## -1.048460e-02 8.040697e-05 -1.279223e-05 7.423443e-06 2.108395e-04
## 1.235831e-05 -1.151088e-07 1.342203e-08 -9.786151e-09 -1.482819e-06
## -4.562025e-05 3.411114e-07 -4.875539e-08 2.841651e-08 7.369382e-06
## 3.739740e-04 -1.424738e-06 3.975327e-07 -1.801662e-07 -2.619119e-04
## -5.974424e-06 3.726223e-08 -6.577606e-09 3.280802e-09 6.340106e-07
## -1.020707e-02 7.829575e-05 -1.236109e-05 7.236109e-06 1.930788e-04
## -1.543032e-06 2.153623e-09 -2.903121e-10 1.500324e-10 6.109753e-09
## 1.744721e-07 -1.010411e-07 1.487890e-09 -4.301793e-10 6.018731e-07
## 6.087068e-08 -6.042358e-10 -1.886722e-09 -3.080407e-11 5.987959e-09
## -6.218570e-07 6.443454e-09 -9.912593e-10 -1.283579e-08 9.193719e-08
## -7.190614e-08 8.399193e-10 -4.535841e-11 5.389332e-11 -5.059743e-07
## -3.658575e-07 3.574915e-09 -6.197751e-10 3.391177e-10 -1.502629e-07
## -3.001037e-08 -1.542040e-11 -5.129840e-11 -2.459479e-12 2.295868e-09
## -4.754007e-08 2.274434e-10 -4.332115e-11 1.868200e-11 2.099079e-08
## -1.124169e-04 9.597573e-07 -1.787641e-07 1.575489e-07 -2.964893e-05
## -4.315669e-08 -1.318156e-10 -1.121372e-10 -1.409411e-11 -4.990925e-08
## -7.895834e-07 6.471339e-09 -9.009602e-10 4.096428e-10 4.468599e-08
## -1.631793e-06 1.827411e-10 1.862549e-10 5.665024e-12 -1.333481e-09
## 2.412921e-06 -7.134720e-06 -7.110455e-08 5.637297e-08 -3.008611e-05
## 5.164326e-06 -4.245534e-08 -3.974910e-07 -3.863068e-09 -1.350963e-06
## -1.059829e-04 8.692692e-07 -2.790784e-07 -3.256117e-06 -5.752260e-06
## -1.285670e-05 3.197720e-08 -1.390769e-07 -4.087233e-08 1.493982e-03
## -2.590908e-05 6.349626e-08 -4.808281e-08 3.239363e-08 -1.533489e-05
## -1.425741e-07 1.262323e-09 -4.233241e-10 3.332600e-11 1.495735e-07
## 2.450263e-06 -1.792318e-08 2.403761e-09 -1.857385e-09 -6.963543e-08
## 1.442059e-04 -4.462257e-07 3.447984e-07 -1.488855e-07 -1.154568e-04
## -2.003394e-04 1.276002e-06 -1.860289e-07 1.488855e-07 -1.898936e-06
## -1.718510e-04 -5.316252e-08 3.804909e-08 1.720440e-07 1.064260e-04
## 2.170644e-03 -1.136000e-07 -3.105157e-07 9.755359e-08 -1.219124e-05
## -1.390778e-05 -2.647988e-06 -2.345355e-07 1.063924e-07 -3.229459e-05
## 1.142873e-04 -6.231211e-07 -6.899279e-06 -9.822213e-09 5.055921e-05
## 9.046515e-05 -6.801336e-07 1.572541e-07 1.502234e-06 -3.899643e-05
## -1.726516e-06 1.089017e-06 -3.947735e-09 -6.297612e-08 2.024029e-03
## 1.044712e-04 -5.095657e-07 -9.850040e-08 -6.635139e-08 -7.633214e-05
## -3.628389e-04 4.890167e-06 -7.032046e-07 3.725968e-07 7.950768e-06
## 1.033828e-06 -9.172890e-09 1.493144e-09 -1.025377e-09 1.978976e-07
## -7.870081e-05 -7.776944e-07 9.306381e-08 -9.599442e-08 1.201458e-05
## 7.555685e-05 -1.348498e-07 -1.344283e-07 1.283345e-08 1.301678e-04
## -6.116661e-05 -4.815620e-07 5.755218e-08 8.671205e-08 1.157309e-04
## 2.257022e-03 -3.071385e-04 8.308783e-06 -5.896391e-05 -3.886636e-02
## 7.199313e-04 2.943031e-06 4.562386e-06 2.194091e-06 -1.755537e-03
## 7.199314e-04 2.943031e-06 4.562386e-06 2.194091e-06 -1.755537e-03
## 5.085231e-03 2.738674e-05 -5.752910e-05 -3.091300e-05 -4.369041e-03
## d -3.082193e-01 -2.020641e-03 -4.052663e-04 -1.414366e-04 -2.431517e+00
## [,6] [,7] [,8] [,9] [,10]
## -1.048460e-02 1.235831e-05 -4.562025e-05 3.739740e-04 -5.974424e-06
## 8.040697e-05 -1.151088e-07 3.411114e-07 -1.424738e-06 3.726223e-08
## -1.279223e-05 1.342203e-08 -4.875539e-08 3.975327e-07 -6.577606e-09
## 7.423443e-06 -9.786151e-09 2.841651e-08 -1.801662e-07 3.280802e-09
## 2.108395e-04 -1.482819e-06 7.369382e-06 -2.619119e-04 6.340106e-07
## 3.009666e-01 7.527102e-06 -2.141337e-05 7.991809e-04 -3.751252e-06
## 7.527102e-06 5.666869e-04 3.209779e-08 -6.330079e-08 3.571425e-09
## -2.141337e-05 3.209779e-08 8.068706e-04 4.375708e-07 -1.644120e-08
## 7.991809e-04 -6.330079e-08 4.375708e-07 8.154803e-01 1.293735e-08
## -3.751252e-06 3.571425e-09 -1.644120e-08 1.293735e-08 3.948809e-04
## -8.990453e-03 7.285399e-06 -2.018448e-05 -4.562371e-04 -3.654175e-06
## -1.613358e-07 8.573457e-11 -7.760877e-10 3.080859e-07 -5.594100e-11
## 5.315091e-07 -2.788568e-11 3.828569e-09 -1.833049e-06 7.862870e-11
## -7.124533e-09 -5.066141e-11 6.913405e-11 9.146134e-09 1.838680e-11
## -2.764332e-07 5.207613e-10 -1.978138e-09 4.409438e-07 -2.842356e-10
## 5.619747e-08 -8.178187e-11 2.787484e-10 6.450662e-08 -5.255105e-12
## 4.402493e-06 2.716184e-10 -8.596090e-10 2.742081e-07 -5.312219e-11
## 4.419438e-10 4.799413e-10 3.909295e-11 3.427410e-09 -2.722563e-13
## -1.782279e-08 4.692699e-11 4.302231e-09 2.068727e-08 -2.092218e-11
## -3.819495e-05 6.820035e-08 -2.848133e-07 -3.203809e-04 -1.012429e-08
## -9.857243e-08 -8.172531e-12 -9.032177e-13 2.025494e-07 -1.157396e-09
## -5.105848e-07 1.829942e-10 -2.516569e-09 8.842014e-07 -1.995331e-10
## 5.868920e-08 -9.962741e-12 -3.003937e-10 1.644142e-07 -4.808818e-12
## -3.877686e-05 -6.637415e-08 1.865045e-07 5.420671e-06 3.506928e-08
## 1.785611e-06 -4.470898e-09 1.981567e-08 -1.111747e-07 1.623018e-09
## -1.130376e-05 8.128656e-08 -2.238067e-07 -1.282022e-04 -6.643348e-09
## -2.138298e-06 2.731278e-08 2.654559e-07 1.147596e-04 4.176476e-08
## 1.428430e-03 -5.612564e-09 2.412665e-07 -8.243374e-05 -1.086992e-08
## 9.359913e-09 -9.150423e-09 -7.453653e-10 -5.326715e-08 2.655058e-12
## 1.533471e-06 -1.958523e-09 -2.524391e-07 -4.680692e-08 9.702083e-10
## 1.888240e-04 -2.698588e-08 1.058543e-06 -1.117623e-02 7.552663e-08
## -6.471768e-05 1.067841e-07 -2.851809e-07 -8.866707e-05 1.600421e-06
## -1.471715e-04 -7.373578e-09 -3.275941e-08 -5.745911e-06 -5.092850e-08
## -5.012560e-05 -7.958430e-08 -1.769486e-07 -1.567767e-04 4.753144e-09
## -1.438868e-05 9.039550e-08 6.768085e-08 -1.344600e-04 -2.434625e-08
## 1.002883e-04 -1.644520e-07 -4.108493e-08 1.932358e-05 8.974287e-09
## 4.943729e-05 -6.970148e-08 1.647441e-07 9.497290e-06 2.792397e-08
## -4.110887e-05 -4.579127e-08 1.899018e-08 1.055461e-04 4.325041e-08
## 4.086243e-04 -1.150133e-07 -1.843548e-07 -1.106863e-04 4.823713e-08
## -1.965383e-04 -9.411508e-06 -1.213808e-06 -1.568175e-04 -1.437105e-07
## 9.388461e-07 -7.213098e-10 -5.200215e-08 3.359290e-07 3.906524e-10
## 5.494384e-07 -5.238045e-08 -4.378520e-07 6.704777e-03 1.067055e-08
## 6.028712e-05 1.225716e-08 3.235285e-07 -1.113359e-04 -1.097967e-06
## 3.425042e-05 1.414281e-07 -1.794619e-07 8.343746e-05 -3.489360e-08
## 4.619034e-02 -2.696429e-05 -7.296186e-05 1.388744e-02 2.490438e-06
## -2.164583e-04 -6.713834e-07 -1.043364e-05 2.861480e-03 -2.029725e-07
## -2.164583e-04 -6.713834e-07 -1.043364e-05 2.861480e-03 -2.029725e-07
## 1.933550e-02 5.993649e-07 -1.490813e-04 -3.019077e-02 -2.342282e-05
## d -6.248768e-01 -1.600701e-04 -1.617446e-03 -2.628047e+00 -3.204051e-05
## [,11] [,12] [,13] [,14] [,15]
## -1.020707e-02 -1.543032e-06 1.744721e-07 6.087068e-08 -6.218570e-07
## 7.829575e-05 2.153623e-09 -1.010411e-07 -6.042358e-10 6.443454e-09
## -1.236109e-05 -2.903121e-10 1.487890e-09 -1.886722e-09 -9.912593e-10
## 7.236109e-06 1.500324e-10 -4.301793e-10 -3.080407e-11 -1.283579e-08
## 1.930788e-04 6.109753e-09 6.018731e-07 5.987959e-09 9.193719e-08
## -8.990453e-03 -1.613358e-07 5.315091e-07 -7.124533e-09 -2.764332e-07
## 7.285399e-06 8.573457e-11 -2.788568e-11 -5.066141e-11 5.207613e-10
## -2.018448e-05 -7.760877e-10 3.828569e-09 6.913405e-11 -1.978138e-09
## -4.562371e-04 3.080859e-07 -1.833049e-06 9.146134e-09 4.409438e-07
## -3.654175e-06 -5.594100e-11 7.862870e-11 1.838680e-11 -2.842356e-10
## 3.367273e-01 -1.739061e-07 1.475869e-06 4.637722e-08 -2.795924e-07
## -1.739061e-07 4.266236e-04 6.786035e-09 1.466041e-10 -1.108583e-09
## 1.475869e-06 6.786035e-09 3.278008e-04 -1.155633e-09 9.254138e-09
## 4.637722e-08 1.466041e-10 -1.155633e-09 3.192208e-04 2.346075e-10
## -2.795924e-07 -1.108583e-09 9.254138e-09 2.346075e-10 7.584674e-04
## -7.010728e-08 -2.419092e-11 -1.445661e-10 1.557981e-11 -2.853727e-10
## -4.853357e-07 -1.609478e-09 1.268210e-08 2.453476e-10 -2.085795e-09
## -2.723623e-08 -2.009474e-10 1.419338e-09 2.684567e-11 -1.716897e-10
## -5.082108e-08 -9.358510e-11 6.898926e-10 1.483190e-11 -1.978511e-10
## -1.299209e-05 -3.118021e-07 9.885951e-07 2.841341e-08 4.402429e-08
## -1.560082e-07 -6.428302e-10 5.996232e-09 8.572461e-11 -5.545445e-10
## 1.256132e-05 -2.417004e-09 1.893326e-08 3.222075e-10 -3.218646e-09
## 4.530669e-08 1.517196e-09 3.968997e-10 4.520876e-12 -4.846273e-11
## -8.524389e-05 -1.287210e-09 5.882713e-06 1.300589e-08 2.297406e-08
## 2.839377e-06 6.643348e-10 -7.879130e-09 -5.036142e-09 6.251558e-10
## -3.488730e-05 1.885660e-09 -3.657953e-08 1.899736e-08 -7.140640e-06
## 4.284427e-05 6.964669e-08 -3.265698e-07 9.346176e-09 1.967903e-08
## 1.315328e-05 6.206410e-08 -4.504425e-07 -2.489708e-09 6.627804e-09
## -2.195166e-08 8.705166e-11 -9.667834e-10 -4.016207e-11 -1.727268e-10
## 2.089913e-06 -7.777394e-10 6.371915e-10 1.650557e-10 -1.063074e-09
## 1.037087e-04 5.612153e-08 -5.447897e-07 9.809128e-09 -5.588796e-08
## -1.064285e-04 4.385122e-08 -1.164076e-07 2.456459e-09 -1.186484e-07
## 8.361692e-03 4.345232e-08 -2.822691e-07 -3.096897e-08 2.868999e-07
## -1.018761e-04 3.241660e-06 3.277188e-07 1.616079e-08 -1.692213e-09
## -3.189194e-05 5.347836e-08 1.952012e-05 4.730867e-08 9.758270e-08
## 1.020567e-04 6.495197e-08 3.799052e-08 8.656827e-07 4.463909e-09
## 7.104153e-05 -5.313580e-08 -4.435104e-08 -8.473817e-09 -4.425609e-06
## 1.162302e-04 -1.805179e-07 5.617243e-07 -8.377932e-09 -8.110249e-08
## 1.922757e-04 8.786742e-10 -2.159896e-07 6.705670e-09 -9.741772e-08
## -2.438801e-04 6.190421e-08 5.767038e-07 4.294756e-08 -4.579977e-08
## 1.022535e-06 8.982863e-11 -1.044585e-09 -4.564347e-11 3.203474e-10
## -3.918064e-05 -2.260282e-07 6.566962e-07 6.489594e-08 -2.539667e-07
## 1.806163e-05 1.412452e-08 -3.508106e-08 1.375701e-08 1.012764e-07
## 1.018976e-03 6.221078e-09 1.007469e-06 1.004202e-08 4.449507e-08
## 4.828186e-02 -3.335379e-05 2.165809e-05 1.014447e-06 1.209976e-05
## -6.871462e-04 -1.411320e-06 -3.621755e-06 1.822038e-06 4.191516e-06
## -6.871463e-04 -1.411320e-06 -3.621755e-06 1.822038e-06 4.191516e-06
## 2.367827e-02 8.519951e-06 8.433820e-05 -2.264135e-05 -2.357830e-05
## d -6.276657e-01 -1.609617e-04 -2.426382e-03 -5.265158e-05 -8.823837e-04
## [,16] [,17] [,18] [,19] [,20]
## -7.190614e-08 -3.658575e-07 -3.001037e-08 -4.754007e-08 -1.124169e-04
## 8.399193e-10 3.574915e-09 -1.542040e-11 2.274434e-10 9.597573e-07
## -4.535841e-11 -6.197751e-10 -5.129840e-11 -4.332115e-11 -1.787641e-07
## 5.389332e-11 3.391177e-10 -2.459479e-12 1.868200e-11 1.575489e-07
## -5.059743e-07 -1.502629e-07 2.295868e-09 2.099079e-08 -2.964893e-05
## 5.619747e-08 4.402493e-06 4.419438e-10 -1.782279e-08 -3.819495e-05
## -8.178187e-11 2.716184e-10 4.799413e-10 4.692699e-11 6.820035e-08
## 2.787484e-10 -8.596090e-10 3.909295e-11 4.302231e-09 -2.848133e-07
## 6.450662e-08 2.742081e-07 3.427410e-09 2.068727e-08 -3.203809e-04
## -5.255105e-12 -5.312219e-11 -2.722563e-13 -2.092218e-11 -1.012429e-08
## -7.010728e-08 -4.853357e-07 -2.723623e-08 -5.082108e-08 -1.299209e-05
## -2.419092e-11 -1.609478e-09 -2.009474e-10 -9.358510e-11 -3.118021e-07
## -1.445661e-10 1.268210e-08 1.419338e-09 6.898926e-10 9.885951e-07
## 1.557981e-11 2.453476e-10 2.684567e-11 1.483190e-11 2.841341e-08
## -2.853727e-10 -2.085795e-09 -1.716897e-10 -1.978511e-10 4.402429e-08
## 7.009843e-04 4.448143e-11 6.206674e-12 -4.499026e-11 6.307436e-08
## 4.448143e-11 7.358724e-04 -2.756374e-10 -1.914776e-10 2.722093e-07
## 6.206674e-12 -2.756374e-10 3.501069e-04 -1.638504e-11 -2.287204e-08
## -4.499026e-11 -1.914776e-10 -1.638504e-11 7.377191e-04 3.355966e-08
## 6.307436e-08 2.722093e-07 -2.287204e-08 3.355966e-08 1.101966e-01
## 1.457283e-10 -1.444218e-09 -1.496290e-10 -6.739506e-11 -1.400585e-07
## 2.704474e-10 -4.576547e-09 -5.381515e-10 -2.252786e-10 -6.013149e-07
## 1.259516e-11 -6.727102e-11 -2.917385e-12 1.469414e-12 2.780821e-08
## -5.092699e-08 8.643127e-08 2.538659e-09 -1.010392e-08 -5.856249e-05
## 1.469609e-10 6.711574e-10 4.317906e-10 1.301324e-10 -8.144748e-07
## -3.505897e-08 1.378341e-07 4.544491e-11 -3.705291e-09 -7.827254e-06
## 5.824708e-07 1.462497e-08 1.063626e-08 -3.624178e-09 5.637313e-06
## -2.823206e-08 -9.806650e-06 -5.346505e-09 -9.279904e-09 5.033815e-05
## -5.785911e-11 3.008580e-10 -3.682899e-11 2.517242e-13 5.784399e-08
## 3.425734e-10 -2.708945e-09 -1.843744e-10 4.711515e-10 2.334976e-07
## 8.152773e-09 -2.900723e-07 2.042343e-08 2.719819e-09 -1.061682e-02
## -3.788315e-08 6.957545e-08 2.487004e-10 -2.490266e-09 2.389199e-05
## -4.913122e-08 9.402409e-08 -7.278938e-09 3.608822e-09 -4.596916e-06
## -2.036865e-07 2.023356e-08 7.255091e-09 1.144163e-08 -1.745202e-05
## -3.140949e-09 -1.186439e-07 1.533777e-08 1.160769e-08 -4.806524e-05
## 1.067658e-08 1.701972e-07 1.249628e-09 -4.240190e-09 -6.273318e-05
## -3.376260e-08 -9.567441e-08 -1.423263e-08 -1.565945e-09 2.545498e-06
## 2.054468e-06 -2.768959e-07 -2.804861e-08 -2.994115e-09 -7.643321e-05
## -1.246306e-08 -2.390910e-06 9.833155e-09 -1.414770e-10 5.521124e-05
## 1.794959e-08 -4.917457e-08 -1.279924e-06 1.928289e-08 1.830923e-05
## -6.579758e-11 1.470772e-10 -3.059411e-11 -3.114329e-09 2.151856e-07
## -4.518821e-08 -7.634662e-08 -4.827067e-09 -3.410721e-08 -1.188022e-02
## 1.070145e-09 -1.809482e-07 -3.656216e-09 1.104674e-09 -5.802650e-06
## 5.349096e-08 -4.423456e-07 -4.180949e-08 3.368711e-10 -4.614490e-05
## -1.144909e-05 4.560059e-05 -6.705592e-06 -3.977057e-06 -2.248659e-03
## -1.630473e-07 -1.111844e-06 -5.396411e-07 -8.855926e-07 -9.376581e-04
## -1.630473e-07 -1.111844e-06 -5.396411e-07 -8.855926e-07 -9.376581e-04
## 6.120716e-06 5.784963e-06 -6.430291e-06 -2.394226e-05 -1.876282e-02
## d -4.039753e-04 -7.420917e-04 -3.642110e-05 -5.600997e-05 -1.259753e+00
## [,21] [,22] [,23] [,24] [,25]
## -4.315669e-08 -7.895834e-07 -1.631793e-06 2.412921e-06 5.164326e-06
## -1.318156e-10 6.471339e-09 1.827411e-10 -7.134720e-06 -4.245534e-08
## -1.121372e-10 -9.009602e-10 1.862549e-10 -7.110455e-08 -3.974910e-07
## -1.409411e-11 4.096428e-10 5.665024e-12 5.637297e-08 -3.863068e-09
## -4.990925e-08 4.468599e-08 -1.333481e-09 -3.008611e-05 -1.350963e-06
## -9.857243e-08 -5.105848e-07 5.868920e-08 -3.877686e-05 1.785611e-06
## -8.172531e-12 1.829942e-10 -9.962741e-12 -6.637415e-08 -4.470898e-09
## -9.032177e-13 -2.516569e-09 -3.003937e-10 1.865045e-07 1.981567e-08
## 2.025494e-07 8.842014e-07 1.644142e-07 5.420671e-06 -1.111747e-07
## -1.157396e-09 -1.995331e-10 -4.808818e-12 3.506928e-08 1.623018e-09
## -1.560082e-07 1.256132e-05 4.530669e-08 -8.524389e-05 2.839377e-06
## -6.428302e-10 -2.417004e-09 1.517196e-09 -1.287210e-09 6.643348e-10
## 5.996232e-09 1.893326e-08 3.968997e-10 5.882713e-06 -7.879130e-09
## 8.572461e-11 3.222075e-10 4.520876e-12 1.300589e-08 -5.036142e-09
## -5.545445e-10 -3.218646e-09 -4.846273e-11 2.297406e-08 6.251558e-10
## 1.457283e-10 2.704474e-10 1.259516e-11 -5.092699e-08 1.469609e-10
## -1.444218e-09 -4.576547e-09 -6.727102e-11 8.643127e-08 6.711574e-10
## -1.496290e-10 -5.381515e-10 -2.917385e-12 2.538659e-09 4.317906e-10
## -6.739506e-11 -2.252786e-10 1.469414e-12 -1.010392e-08 1.301324e-10
## -1.400585e-07 -6.013149e-07 2.780821e-08 -5.856249e-05 -8.144748e-07
## 4.371300e-04 -2.499489e-09 -2.967340e-11 -4.380030e-08 1.499192e-10
## -2.499489e-09 4.398353e-04 -8.471994e-11 -2.723448e-07 1.916007e-09
## -2.967340e-11 -8.471994e-11 4.869405e-04 8.991525e-08 -5.724951e-10
## -4.380030e-08 -2.723448e-07 8.991525e-08 2.338156e-01 7.376277e-07
## 1.499192e-10 1.916007e-09 -5.724951e-10 7.376277e-07 4.997102e-04
## 1.119052e-07 -2.997126e-09 9.151215e-08 -7.130569e-05 2.761427e-07
## 5.180879e-08 8.732144e-08 1.158235e-07 -8.562892e-05 1.273263e-06
## 4.002767e-08 1.212859e-07 7.441252e-08 -4.542706e-05 3.882098e-07
## 1.564309e-10 3.347789e-11 1.554233e-10 -7.739377e-08 3.530020e-09
## 8.133167e-10 -7.415015e-10 -7.803963e-10 6.300487e-07 -6.158563e-09
## 3.073715e-08 -5.071929e-07 4.915397e-08 -6.786429e-05 1.215871e-06
## 5.204937e-07 -5.686562e-08 9.517065e-08 -4.485866e-05 1.085941e-06
## 2.916176e-08 -8.385712e-07 -1.669098e-08 -3.452247e-05 9.074281e-07
## -1.375813e-07 -5.519414e-08 1.226598e-05 1.719235e-04 -3.736003e-06
## -7.014377e-08 -9.738179e-09 -1.689676e-08 2.391633e-03 -9.585003e-08
## 2.909270e-08 -3.733214e-08 -3.122089e-08 5.156948e-05 4.265170e-04
## 9.589976e-08 4.021509e-08 1.869053e-08 2.231088e-05 9.402151e-07
## -8.707346e-08 -7.547534e-07 3.732639e-08 -2.189275e-05 -3.467485e-07
## 1.280466e-07 2.284189e-07 -2.057037e-08 7.410912e-05 -9.622106e-07
## -1.545261e-08 -1.796102e-07 4.207362e-09 4.830347e-05 -1.307917e-07
## -3.765369e-11 -7.660971e-11 -1.612689e-11 7.685923e-08 2.147797e-09
## 7.929069e-08 -4.624256e-07 7.278180e-08 1.026267e-04 -3.023369e-06
## -3.386227e-06 -3.387419e-08 4.019895e-08 1.855503e-05 -1.071513e-07
## 1.287187e-07 1.149054e-05 7.924567e-08 1.985897e-05 1.634432e-07
## -8.976204e-06 -5.375594e-05 1.339549e-06 1.373922e-03 5.790584e-06
## 3.570917e-07 -3.256718e-06 1.621304e-07 -6.210512e-04 2.160833e-04
## 3.570917e-07 -3.256718e-06 1.621304e-07 -6.210512e-04 2.160833e-04
## 2.706071e-05 7.068462e-05 1.563087e-06 -6.922986e-04 -2.796928e-03
## d -2.546047e-04 -1.078654e-03 -1.166586e-04 -6.392809e-01 -8.803847e-03
## [,26] [,27] [,28] [,29] [,30]
## -1.059829e-04 -1.285670e-05 -2.590908e-05 -1.425741e-07 2.450263e-06
## 8.692692e-07 3.197720e-08 6.349626e-08 1.262323e-09 -1.792318e-08
## -2.790784e-07 -1.390769e-07 -4.808281e-08 -4.233241e-10 2.403761e-09
## -3.256117e-06 -4.087233e-08 3.239363e-08 3.332600e-11 -1.857385e-09
## -5.752260e-06 1.493982e-03 -1.533489e-05 1.495735e-07 -6.963543e-08
## -1.130376e-05 -2.138298e-06 1.428430e-03 9.359913e-09 1.533471e-06
## 8.128656e-08 2.731278e-08 -5.612564e-09 -9.150423e-09 -1.958523e-09
## -2.238067e-07 2.654559e-07 2.412665e-07 -7.453653e-10 -2.524391e-07
## -1.282022e-04 1.147596e-04 -8.243374e-05 -5.326715e-08 -4.680692e-08
## -6.643348e-09 4.176476e-08 -1.086992e-08 2.655058e-12 9.702083e-10
## -3.488730e-05 4.284427e-05 1.315328e-05 -2.195166e-08 2.089913e-06
## 1.885660e-09 6.964669e-08 6.206410e-08 8.705166e-11 -7.777394e-10
## -3.657953e-08 -3.265698e-07 -4.504425e-07 -9.667834e-10 6.371915e-10
## 1.899736e-08 9.346176e-09 -2.489708e-09 -4.016207e-11 1.650557e-10
## -7.140640e-06 1.967903e-08 6.627804e-09 -1.727268e-10 -1.063074e-09
## -3.505897e-08 5.824708e-07 -2.823206e-08 -5.785911e-11 3.425734e-10
## 1.378341e-07 1.462497e-08 -9.806650e-06 3.008580e-10 -2.708945e-09
## 4.544491e-11 1.063626e-08 -5.346505e-09 -3.682899e-11 -1.843744e-10
## -3.705291e-09 -3.624178e-09 -9.279904e-09 2.517242e-13 4.711515e-10
## -7.827254e-06 5.637313e-06 5.033815e-05 5.784399e-08 2.334976e-07
## 1.119052e-07 5.180879e-08 4.002767e-08 1.564309e-10 8.133167e-10
## -2.997126e-09 8.732144e-08 1.212859e-07 3.347789e-11 -7.415015e-10
## 9.151215e-08 1.158235e-07 7.441252e-08 1.554233e-10 -7.803963e-10
## -7.130569e-05 -8.562892e-05 -4.542706e-05 -7.739377e-08 6.300487e-07
## 2.761427e-07 1.273263e-06 3.882098e-07 3.530020e-09 -6.158563e-09
## 1.470224e-01 -2.736697e-05 -3.339552e-05 -9.009782e-08 3.182207e-07
## -2.736697e-05 1.358293e+00 -4.117275e-05 -1.554863e-07 5.813950e-07
## -3.339552e-05 -4.117275e-05 1.671927e+00 -5.689658e-08 1.204765e-06
## -9.009782e-08 -1.554863e-07 -5.689658e-08 7.956294e-04 8.033030e-10
## 3.182207e-07 5.813950e-07 1.204765e-06 8.033030e-10 1.628934e-05
## 4.008510e-06 1.016749e-05 8.371287e-06 -1.771898e-07 -4.370590e-07
## -4.367812e-05 -5.830189e-05 -5.084189e-05 -1.096259e-07 4.635592e-07
## 5.737395e-05 -7.737827e-05 -7.473019e-05 4.781551e-08 9.028057e-07
## -7.897256e-05 5.226408e-05 6.653994e-05 7.267183e-08 -6.662087e-07
## -2.148782e-05 1.549085e-05 4.572824e-05 3.349571e-08 -4.921022e-07
## 2.073853e-05 2.642033e-06 -2.012045e-05 1.035229e-08 -4.543993e-08
## -3.096115e-02 -4.289334e-06 -2.095196e-05 6.866574e-08 3.295940e-08
## -4.633749e-05 -4.690785e-03 4.144324e-05 3.660438e-08 -3.937636e-07
## 5.758625e-05 5.987719e-05 4.998875e-03 7.545791e-08 -4.282192e-07
## 2.207420e-06 -2.605378e-05 2.826307e-05 -5.049141e-05 -3.647371e-07
## 1.002223e-07 -1.236488e-08 3.704119e-08 2.958895e-11 1.126303e-07
## 1.362103e-04 4.602808e-05 -8.328535e-06 -1.845418e-08 -3.907396e-07
## -4.139481e-05 2.616724e-05 -5.405265e-05 -2.184794e-08 -2.916434e-07
## -1.760852e-05 1.022853e-05 -4.281409e-05 -2.246206e-08 1.813971e-07
## -1.028977e-02 -1.872117e-02 -3.001015e-03 -5.233021e-06 -9.029216e-05
## -2.995065e-03 -2.791267e-04 -8.469429e-04 1.411897e-07 9.008180e-06
## -2.995065e-03 -2.791267e-04 -8.469429e-04 1.411897e-07 9.008180e-06
## 3.703224e-02 9.963481e-03 -1.366001e-02 2.203648e-05 1.346345e-03
## d -4.565047e-01 -1.366577e+00 -1.443090e+00 -2.452335e-04 -2.632642e-03
## [,31] [,32] [,33] [,34] [,35]
## 1.442059e-04 -2.003394e-04 -1.718510e-04 2.170644e-03 -1.390778e-05
## -4.462257e-07 1.276002e-06 -5.316252e-08 -1.136000e-07 -2.647988e-06
## 3.447984e-07 -1.860289e-07 3.804909e-08 -3.105157e-07 -2.345355e-07
## -1.488855e-07 1.488855e-07 1.720440e-07 9.755359e-08 1.063924e-07
## -1.154568e-04 -1.898936e-06 1.064260e-04 -1.219124e-05 -3.229459e-05
## 1.888240e-04 -6.471768e-05 -1.471715e-04 -5.012560e-05 -1.438868e-05
## -2.698588e-08 1.067841e-07 -7.373578e-09 -7.958430e-08 9.039550e-08
## 1.058543e-06 -2.851809e-07 -3.275941e-08 -1.769486e-07 6.768085e-08
## -1.117623e-02 -8.866707e-05 -5.745911e-06 -1.567767e-04 -1.344600e-04
## 7.552663e-08 1.600421e-06 -5.092850e-08 4.753144e-09 -2.434625e-08
## 1.037087e-04 -1.064285e-04 8.361692e-03 -1.018761e-04 -3.189194e-05
## 5.612153e-08 4.385122e-08 4.345232e-08 3.241660e-06 5.347836e-08
## -5.447897e-07 -1.164076e-07 -2.822691e-07 3.277188e-07 1.952012e-05
## 9.809128e-09 2.456459e-09 -3.096897e-08 1.616079e-08 4.730867e-08
## -5.588796e-08 -1.186484e-07 2.868999e-07 -1.692213e-09 9.758270e-08
## 8.152773e-09 -3.788315e-08 -4.913122e-08 -2.036865e-07 -3.140949e-09
## -2.900723e-07 6.957545e-08 9.402409e-08 2.023356e-08 -1.186439e-07
## 2.042343e-08 2.487004e-10 -7.278938e-09 7.255091e-09 1.533777e-08
## 2.719819e-09 -2.490266e-09 3.608822e-09 1.144163e-08 1.160769e-08
## -1.061682e-02 2.389199e-05 -4.596916e-06 -1.745202e-05 -4.806524e-05
## 3.073715e-08 5.204937e-07 2.916176e-08 -1.375813e-07 -7.014377e-08
## -5.071929e-07 -5.686562e-08 -8.385712e-07 -5.519414e-08 -9.738179e-09
## 4.915397e-08 9.517065e-08 -1.669098e-08 1.226598e-05 -1.689676e-08
## -6.786429e-05 -4.485866e-05 -3.452247e-05 1.719235e-04 2.391633e-03
## 1.215871e-06 1.085941e-06 9.074281e-07 -3.736003e-06 -9.585003e-08
## 4.008510e-06 -4.367812e-05 5.737395e-05 -7.897256e-05 -2.148782e-05
## 1.016749e-05 -5.830189e-05 -7.737827e-05 5.226408e-05 1.549085e-05
## 8.371287e-06 -5.084189e-05 -7.473019e-05 6.653994e-05 4.572824e-05
## -1.771898e-07 -1.096259e-07 4.781551e-08 7.267183e-08 3.349571e-08
## -4.370590e-07 4.635592e-07 9.028057e-07 -6.662087e-07 -4.921022e-07
## 8.794948e-01 2.149916e-06 8.098385e-05 1.057290e-04 7.042952e-05
## 2.149916e-06 9.114423e-02 1.519498e-05 7.319854e-05 6.009655e-05
## 8.098385e-05 1.519498e-05 4.711383e-01 2.482397e-05 1.021419e-04
## 1.057290e-04 7.319854e-05 2.482397e-05 8.724744e-01 -2.680573e-04
## 7.042952e-05 6.009655e-05 1.021419e-04 -2.680573e-04 1.559592e-01
## 4.141658e-06 1.082377e-05 4.889460e-05 -2.067287e-05 -7.938970e-05
## -6.457024e-05 1.743429e-05 -2.302713e-05 -9.473348e-05 -6.523228e-05
## -3.941327e-05 5.797262e-06 1.215868e-07 1.263696e-04 -4.898164e-05
## 1.727148e-04 -2.335329e-06 4.418983e-05 -2.042897e-04 -1.346788e-04
## 2.532689e-05 -2.917503e-05 -1.384580e-05 -2.008798e-04 -2.317168e-04
## 5.266079e-08 -5.075766e-09 -9.741729e-08 4.588320e-07 3.053105e-07
## 6.896682e-03 -5.536626e-05 -9.169954e-05 -1.076787e-05 2.059903e-07
## -2.710148e-05 -2.805423e-02 1.492911e-05 -5.870660e-05 -3.726048e-05
## -6.393389e-06 -3.762828e-05 -2.491598e-03 1.877426e-04 -2.013059e-05
## -1.135529e-02 -1.621217e-02 4.825363e-03 2.049162e-02 2.804752e-02
## 2.502683e-03 8.834319e-04 9.443769e-04 -9.754980e-04 7.703872e-04
## 2.502683e-03 8.834319e-04 9.443769e-04 -9.754981e-04 7.703872e-04
## 1.055105e-03 -6.011211e-03 5.571056e-03 3.258315e-03 9.493410e-03
## d -2.237778e+00 -2.792301e-01 -1.301326e+00 -1.418846e+00 -6.990457e-01
## [,36] [,37] [,38] [,39] [,40]
## 1.142873e-04 9.046515e-05 -1.726516e-06 1.044712e-04 -3.628389e-04
## -6.231211e-07 -6.801336e-07 1.089017e-06 -5.095657e-07 4.890167e-06
## -6.899279e-06 1.572541e-07 -3.947735e-09 -9.850040e-08 -7.032046e-07
## -9.822213e-09 1.502234e-06 -6.297612e-08 -6.635139e-08 3.725968e-07
## 5.055921e-05 -3.899643e-05 2.024029e-03 -7.633214e-05 7.950768e-06
## 1.002883e-04 4.943729e-05 -4.110887e-05 4.086243e-04 -1.965383e-04
## -1.644520e-07 -6.970148e-08 -4.579127e-08 -1.150133e-07 -9.411508e-06
## -4.108493e-08 1.647441e-07 1.899018e-08 -1.843548e-07 -1.213808e-06
## 1.932358e-05 9.497290e-06 1.055461e-04 -1.106863e-04 -1.568175e-04
## 8.974287e-09 2.792397e-08 4.325041e-08 4.823713e-08 -1.437105e-07
## 1.020567e-04 7.104153e-05 1.162302e-04 1.922757e-04 -2.438801e-04
## 6.495197e-08 -5.313580e-08 -1.805179e-07 8.786742e-10 6.190421e-08
## 3.799052e-08 -4.435104e-08 5.617243e-07 -2.159896e-07 5.767038e-07
## 8.656827e-07 -8.473817e-09 -8.377932e-09 6.705670e-09 4.294756e-08
## 4.463909e-09 -4.425609e-06 -8.110249e-08 -9.741772e-08 -4.579977e-08
## 1.067658e-08 -3.376260e-08 2.054468e-06 -1.246306e-08 1.794959e-08
## 1.701972e-07 -9.567441e-08 -2.768959e-07 -2.390910e-06 -4.917457e-08
## 1.249628e-09 -1.423263e-08 -2.804861e-08 9.833155e-09 -1.279924e-06
## -4.240190e-09 -1.565945e-09 -2.994115e-09 -1.414770e-10 1.928289e-08
## -6.273318e-05 2.545498e-06 -7.643321e-05 5.521124e-05 1.830923e-05
## 2.909270e-08 9.589976e-08 -8.707346e-08 1.280466e-07 -1.545261e-08
## -3.733214e-08 4.021509e-08 -7.547534e-07 2.284189e-07 -1.796102e-07
## -3.122089e-08 1.869053e-08 3.732639e-08 -2.057037e-08 4.207362e-09
## 5.156948e-05 2.231088e-05 -2.189275e-05 7.410912e-05 4.830347e-05
## 4.265170e-04 9.402151e-07 -3.467485e-07 -9.622106e-07 -1.307917e-07
## 2.073853e-05 -3.096115e-02 -4.633749e-05 5.758625e-05 2.207420e-06
## 2.642033e-06 -4.289334e-06 -4.690785e-03 5.987719e-05 -2.605378e-05
## -2.012045e-05 -2.095196e-05 4.144324e-05 4.998875e-03 2.826307e-05
## 1.035229e-08 6.866574e-08 3.660438e-08 7.545791e-08 -5.049141e-05
## -4.543993e-08 3.295940e-08 -3.937636e-07 -4.282192e-07 -3.647371e-07
## 4.141658e-06 -6.457024e-05 -3.941327e-05 1.727148e-04 2.532689e-05
## 1.082377e-05 1.743429e-05 5.797262e-06 -2.335329e-06 -2.917503e-05
## 4.889460e-05 -2.302713e-05 1.215868e-07 4.418983e-05 -1.384580e-05
## -2.067287e-05 -9.473348e-05 1.263696e-04 -2.042897e-04 -2.008798e-04
## -7.938970e-05 -6.523228e-05 -4.898164e-05 -1.346788e-04 -2.317168e-04
## 1.009182e+00 -6.298454e-05 7.520988e-05 -1.240703e-04 -7.520446e-05
## -6.298454e-05 8.620389e-01 4.364974e-05 -5.402793e-05 -1.308903e-04
## 7.520988e-05 4.364974e-05 1.306747e+00 3.109113e-04 3.403742e-05
## -1.240703e-04 -5.402793e-05 3.109113e-04 1.654121e+00 -1.671879e-04
## -7.520446e-05 -1.308903e-04 3.403742e-05 -1.671879e-04 8.202010e-01
## 1.992967e-07 3.053799e-07 -5.726571e-07 5.890778e-07 3.886130e-07
## -6.678670e-05 -1.921869e-05 -2.401028e-04 -2.234605e-04 -1.496755e-04
## -9.479824e-05 -7.476932e-05 3.512099e-05 -3.899937e-05 -7.133745e-05
## -1.059712e-04 -1.370657e-04 -1.327184e-05 -2.469326e-05 -7.022160e-05
## -7.395783e-03 4.780024e-03 -3.013235e-02 -1.447738e-02 1.949866e-02
## -6.886543e-04 -1.345697e-04 2.533807e-04 2.060231e-03 -2.886398e-05
## -6.886543e-04 -1.345697e-04 2.533807e-04 2.060231e-03 -2.886398e-05
## 6.854026e-02 3.302694e-03 1.143910e-02 8.089496e-02 6.330225e-02
## d -1.138453e+00 -1.016551e+00 -2.406411e+00 -1.920326e+00 -1.059786e+00
## [,41] [,42] [,43] [,44] [,45]
## 1.033828e-06 -7.870081e-05 7.555685e-05 -6.116661e-05 2.257022e-03
## -9.172890e-09 -7.776944e-07 -1.348498e-07 -4.815620e-07 -3.071385e-04
## 1.493144e-09 9.306381e-08 -1.344283e-07 5.755218e-08 8.308783e-06
## -1.025377e-09 -9.599442e-08 1.283345e-08 8.671205e-08 -5.896391e-05
## 1.978976e-07 1.201458e-05 1.301678e-04 1.157309e-04 -3.886636e-02
## 9.388461e-07 5.494384e-07 6.028712e-05 3.425042e-05 4.619034e-02
## -7.213098e-10 -5.238045e-08 1.225716e-08 1.414281e-07 -2.696429e-05
## -5.200215e-08 -4.378520e-07 3.235285e-07 -1.794619e-07 -7.296186e-05
## 3.359290e-07 6.704777e-03 -1.113359e-04 8.343746e-05 1.388744e-02
## 3.906524e-10 1.067055e-08 -1.097967e-06 -3.489360e-08 2.490438e-06
## 1.022535e-06 -3.918064e-05 1.806163e-05 1.018976e-03 4.828186e-02
## 8.982863e-11 -2.260282e-07 1.412452e-08 6.221078e-09 -3.335379e-05
## -1.044585e-09 6.566962e-07 -3.508106e-08 1.007469e-06 2.165809e-05
## -4.564347e-11 6.489594e-08 1.375701e-08 1.004202e-08 1.014447e-06
## 3.203474e-10 -2.539667e-07 1.012764e-07 4.449507e-08 1.209976e-05
## -6.579758e-11 -4.518821e-08 1.070145e-09 5.349096e-08 -1.144909e-05
## 1.470772e-10 -7.634662e-08 -1.809482e-07 -4.423456e-07 4.560059e-05
## -3.059411e-11 -4.827067e-09 -3.656216e-09 -4.180949e-08 -6.705592e-06
## -3.114329e-09 -3.410721e-08 1.104674e-09 3.368711e-10 -3.977057e-06
## 2.151856e-07 -1.188022e-02 -5.802650e-06 -4.614490e-05 -2.248659e-03
## -3.765369e-11 7.929069e-08 -3.386227e-06 1.287187e-07 -8.976204e-06
## -7.660971e-11 -4.624256e-07 -3.387419e-08 1.149054e-05 -5.375594e-05
## -1.612689e-11 7.278180e-08 4.019895e-08 7.924567e-08 1.339549e-06
## 7.685923e-08 1.026267e-04 1.855503e-05 1.985897e-05 1.373922e-03
## 2.147797e-09 -3.023369e-06 -1.071513e-07 1.634432e-07 5.790584e-06
## 1.002223e-07 1.362103e-04 -4.139481e-05 -1.760852e-05 -1.028977e-02
## -1.236488e-08 4.602808e-05 2.616724e-05 1.022853e-05 -1.872117e-02
## 3.704119e-08 -8.328535e-06 -5.405265e-05 -4.281409e-05 -3.001015e-03
## 2.958895e-11 -1.845418e-08 -2.184794e-08 -2.246206e-08 -5.233021e-06
## 1.126303e-07 -3.907396e-07 -2.916434e-07 1.813971e-07 -9.029216e-05
## 5.266079e-08 6.896682e-03 -2.710148e-05 -6.393389e-06 -1.135529e-02
## -5.075766e-09 -5.536626e-05 -2.805423e-02 -3.762828e-05 -1.621217e-02
## -9.741729e-08 -9.169954e-05 1.492911e-05 -2.491598e-03 4.825363e-03
## 4.588320e-07 -1.076787e-05 -5.870660e-05 1.877426e-04 2.049162e-02
## 3.053105e-07 2.059903e-07 -3.726048e-05 -2.013059e-05 2.804752e-02
## 1.992967e-07 -6.678670e-05 -9.479824e-05 -1.059712e-04 -7.395783e-03
## 3.053799e-07 -1.921869e-05 -7.476932e-05 -1.370657e-04 4.780024e-03
## -5.726571e-07 -2.401028e-04 3.512099e-05 -1.327184e-05 -3.013235e-02
## 5.890778e-07 -2.234605e-04 -3.899937e-05 -2.469326e-05 -1.447738e-02
## 3.886130e-07 -1.496755e-04 -7.133745e-05 -7.022160e-05 1.949866e-02
## 4.679136e-04 -1.784104e-07 2.222265e-07 5.411390e-08 -5.032380e-05
## -1.784104e-07 2.284735e+00 6.117901e-05 -7.190204e-05 1.640786e-02
## 2.222265e-07 6.117901e-05 1.038815e+00 -5.538334e-05 3.433426e-03
## 5.411390e-08 -7.190204e-05 -5.538334e-05 3.691314e-01 9.310452e-05
## -5.032380e-05 1.640786e-02 3.433426e-03 9.310452e-05 7.852636e+01
## -1.322123e-05 -1.974175e-03 1.956707e-03 -5.626559e-04 7.194970e-02
## -1.322123e-05 -1.974175e-03 1.956707e-03 -5.626559e-04 7.194970e-02
## 3.977290e-04 -1.465292e-01 2.604409e-02 1.175368e-03 1.037284e+00
## d -1.273815e-03 -3.429783e+00 -1.031776e+00 -1.419134e+00 -1.644841e+02
## [,46] [,47] [,48] [,49]
## 7.199313e-04 7.199314e-04 5.085231e-03 -3.082193e-01
## 2.943031e-06 2.943031e-06 2.738674e-05 -2.020641e-03
## 4.562386e-06 4.562386e-06 -5.752910e-05 -4.052663e-04
## 2.194091e-06 2.194091e-06 -3.091300e-05 -1.414366e-04
## -1.755537e-03 -1.755537e-03 -4.369041e-03 -2.431517e+00
## -2.164583e-04 -2.164583e-04 1.933550e-02 -6.248768e-01
## -6.713834e-07 -6.713834e-07 5.993649e-07 -1.600701e-04
## -1.043364e-05 -1.043364e-05 -1.490813e-04 -1.617446e-03
## 2.861480e-03 2.861480e-03 -3.019077e-02 -2.628047e+00
## -2.029725e-07 -2.029725e-07 -2.342282e-05 -3.204051e-05
## -6.871462e-04 -6.871463e-04 2.367827e-02 -6.276657e-01
## -1.411320e-06 -1.411320e-06 8.519951e-06 -1.609617e-04
## -3.621755e-06 -3.621755e-06 8.433820e-05 -2.426382e-03
## 1.822038e-06 1.822038e-06 -2.264135e-05 -5.265158e-05
## 4.191516e-06 4.191516e-06 -2.357830e-05 -8.823837e-04
## -1.630473e-07 -1.630473e-07 6.120716e-06 -4.039753e-04
## -1.111844e-06 -1.111844e-06 5.784963e-06 -7.420917e-04
## -5.396411e-07 -5.396411e-07 -6.430291e-06 -3.642110e-05
## -8.855926e-07 -8.855926e-07 -2.394226e-05 -5.600997e-05
## -9.376581e-04 -9.376581e-04 -1.876282e-02 -1.259753e+00
## 3.570917e-07 3.570917e-07 2.706071e-05 -2.546047e-04
## -3.256718e-06 -3.256718e-06 7.068462e-05 -1.078654e-03
## 1.621304e-07 1.621304e-07 1.563087e-06 -1.166586e-04
## -6.210512e-04 -6.210512e-04 -6.922986e-04 -6.392809e-01
## 2.160833e-04 2.160833e-04 -2.796928e-03 -8.803847e-03
## -2.995065e-03 -2.995065e-03 3.703224e-02 -4.565047e-01
## -2.791267e-04 -2.791267e-04 9.963481e-03 -1.366577e+00
## -8.469429e-04 -8.469429e-04 -1.366001e-02 -1.443090e+00
## 1.411897e-07 1.411897e-07 2.203648e-05 -2.452335e-04
## 9.008180e-06 9.008180e-06 1.346345e-03 -2.632642e-03
## 2.502683e-03 2.502683e-03 1.055105e-03 -2.237778e+00
## 8.834319e-04 8.834319e-04 -6.011211e-03 -2.792301e-01
## 9.443769e-04 9.443769e-04 5.571056e-03 -1.301326e+00
## -9.754980e-04 -9.754981e-04 3.258315e-03 -1.418846e+00
## 7.703872e-04 7.703872e-04 9.493410e-03 -6.990457e-01
## -6.886543e-04 -6.886543e-04 6.854026e-02 -1.138453e+00
## -1.345697e-04 -1.345697e-04 3.302694e-03 -1.016551e+00
## 2.533807e-04 2.533807e-04 1.143910e-02 -2.406411e+00
## 2.060231e-03 2.060231e-03 8.089496e-02 -1.920326e+00
## -2.886398e-05 -2.886398e-05 6.330225e-02 -1.059786e+00
## -1.322123e-05 -1.322123e-05 3.977290e-04 -1.273815e-03
## -1.974175e-03 -1.974175e-03 -1.465292e-01 -3.429783e+00
## 1.956707e-03 1.956707e-03 2.604409e-02 -1.031776e+00
## -5.626559e-04 -5.626559e-04 1.175368e-03 -1.419134e+00
## 7.194970e-02 7.194970e-02 1.037284e+00 -1.644841e+02
## 2.150837e+01 2.150837e+01 4.453800e+00 -4.804158e+01
## 2.150837e+01 2.150837e+01 4.453800e+00 -4.804158e+01
## 4.453800e+00 4.453800e+00 8.352199e+02 -9.571858e+02
## d -4.804158e+01 -4.804158e+01 -9.571858e+02 5.464350e+04
##
## Model rank = 4080 / 4080
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time):VowelDRESS 9.00 1.62 0.99 0.35
## s(speaker_scaled_time):VowelFLEECE 9.00 1.01 0.99 0.33
## s(speaker_scaled_time):VowelFOOT 9.00 1.00 0.99 0.33
## s(speaker_scaled_time):VowelGOOSE 9.00 1.00 0.99 0.32
## s(speaker_scaled_time):VowelKIT 9.00 5.86 0.99 0.30
## s(speaker_scaled_time):VowelLOT 9.00 2.25 0.99 0.36
## s(speaker_scaled_time):VowelNURSE 9.00 1.00 0.99 0.30
## s(speaker_scaled_time):VowelSTART 9.00 1.00 0.99 0.32
## s(speaker_scaled_time):VowelSTRUT 9.00 6.26 0.99 0.28
## s(speaker_scaled_time):VowelTHOUGHT 9.00 1.00 0.99 0.30
## s(speaker_scaled_time):VowelTRAP 9.00 2.26 0.99 0.38
## s(speaker_scaled_art_rate):VowelDRESS 9.00 1.00 1.00 0.34
## s(speaker_scaled_art_rate):VowelFLEECE 9.00 1.01 1.00 0.45
## s(speaker_scaled_art_rate):VowelFOOT 9.00 1.00 1.00 0.46
## s(speaker_scaled_art_rate):VowelGOOSE 9.00 1.00 1.00 0.45
## s(speaker_scaled_art_rate):VowelKIT 9.00 1.00 1.00 0.41
## s(speaker_scaled_art_rate):VowelLOT 9.00 1.00 1.00 0.39
## s(speaker_scaled_art_rate):VowelNURSE 9.00 1.00 1.00 0.38
## s(speaker_scaled_art_rate):VowelSTART 9.00 1.00 1.00 0.48
## s(speaker_scaled_art_rate):VowelSTRUT 9.00 3.52 1.00 0.42
## s(speaker_scaled_art_rate):VowelTHOUGHT 9.00 1.00 1.00 0.51
## s(speaker_scaled_art_rate):VowelTRAP 9.00 1.00 1.00 0.42
## s(speaker_scaled_amp_max):VowelDRESS 9.00 1.00 0.98 0.17
## s(speaker_scaled_amp_max):VowelFLEECE 9.00 2.28 0.98 0.17
## s(speaker_scaled_amp_max):VowelFOOT 9.00 1.02 0.98 0.17
## s(speaker_scaled_amp_max):VowelGOOSE 9.00 1.91 0.98 0.11
## s(speaker_scaled_amp_max):VowelKIT 9.00 3.73 0.98 0.13
## s(speaker_scaled_amp_max):VowelLOT 9.00 3.89 0.98 0.16
## s(speaker_scaled_amp_max):VowelNURSE 9.00 1.00 0.98 0.16
## s(speaker_scaled_amp_max):VowelSTART 9.00 1.01 0.98 0.13
## s(speaker_scaled_amp_max):VowelSTRUT 9.00 5.48 0.98 0.14
## s(speaker_scaled_amp_max):VowelTHOUGHT 9.00 1.56 0.98 0.14
## s(speaker_scaled_amp_max):VowelTRAP 9.00 3.60 0.98 0.13
## s(speaker_scaled_pitch):VowelDRESS 9.00 3.84 1.00 0.46
## s(speaker_scaled_pitch):VowelFLEECE 9.00 2.40 1.00 0.44
## s(speaker_scaled_pitch):VowelFOOT 9.00 3.28 1.00 0.37
## s(speaker_scaled_pitch):VowelGOOSE 9.00 3.03 1.00 0.45
## s(speaker_scaled_pitch):VowelKIT 9.00 5.81 1.00 0.43
## s(speaker_scaled_pitch):VowelLOT 9.00 4.84 1.00 0.43
## s(speaker_scaled_pitch):VowelNURSE 9.00 3.12 1.00 0.42
## s(speaker_scaled_pitch):VowelSTART 9.00 1.00 1.00 0.44
## s(speaker_scaled_pitch):VowelSTRUT 9.00 7.86 1.00 0.49
## s(speaker_scaled_pitch):VowelTHOUGHT 9.00 3.06 1.00 0.46
## s(speaker_scaled_pitch):VowelTRAP 9.00 3.84 1.00 0.42
## s(speaker_scaled_time,Speaker) 1080.00 425.05 0.99 0.33
## s(Speaker) 216.00 96.08 NA NA
## s(Vowel,Speaker) 2376.00 1914.37 NA NA
No great gain is achieved here. The model output is remarkably similar for our main variable on interest (amplitude). Whatever improvement we might be getting is undermined by the increased computational cost.
The QQ-plots and residual histograms for the previous two structures show that
our residuals do not exactly follow a normal distribution. They have heavy
tails. One way to handle heavy tailed distributions of residuals is to fit a
model which does not assume the normality of the residuals. One such distribution
is the scaled-t distribution (Sóskuthy 2021, 20).
However, fitting with a scaled t distribution (using family = scat(link='identity'))
increased the computational resource demands of the model.
gamm_fit_scat <- bam(
F1_50 ~
participant_gender +
Vowel +
s(speaker_scaled_time, by=Vowel) +
s(speaker_scaled_art_rate, by=Vowel) +
s(speaker_scaled_amp_max, by=Vowel) +
s(speaker_scaled_pitch, by=Vowel) +
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="fREML",
discrete = TRUE,
family = scat(link='identity'),
nthreads = 8 # Change this depending on number of cores available.
)
write_rds(gamm_fit_scat, here('models', 'gamm_fit_scat.rds'))
# calculate summary and save.
gamm_fit_scat_summary <- summary(gamm_fit_scat)
write_rds(gamm_fit_scat_summary, here('models', 'gamm_fit_scat_summary.rds'))
We look at the model summary.
gamm_fit_scat <- read_rds(here('models', 'gamm_fit_scat.rds'))
gamm_fit_scat_summary <- read_rds(here('models', 'gamm_fit_scat_summary.rds'))
gamm_fit_scat_summary
##
## Family: Scaled t(3.428,48.361)
## Link function: identity
##
## Formula:
## F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(Speaker,
## bs = "re") + s(Speaker, Vowel, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 432.520 4.109 105.251 < 2e-16 ***
## participant_genderM -55.497 5.610 -9.893 < 2e-16 ***
## VowelFLEECE -13.632 4.002 -3.406 0.000658 ***
## VowelFOOT 65.225 4.368 14.932 < 2e-16 ***
## VowelGOOSE -21.258 4.140 -5.135 2.83e-07 ***
## VowelKIT 83.824 4.022 20.841 < 2e-16 ***
## VowelLOT 163.255 4.093 39.884 < 2e-16 ***
## VowelNURSE 16.299 4.080 3.995 6.47e-05 ***
## VowelSTART 363.020 4.215 86.123 < 2e-16 ***
## VowelSTRUT 255.864 4.150 61.661 < 2e-16 ***
## VowelTHOUGHT 11.858 4.030 2.943 0.003256 **
## VowelTRAP 146.833 4.107 35.756 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(speaker_scaled_time):VowelDRESS 2.151 2.685 4.171 0.009177 **
## s(speaker_scaled_time):VowelFLEECE 1.077 1.150 4.919 0.025880 *
## s(speaker_scaled_time):VowelFOOT 1.003 1.006 3.693 0.054163 .
## s(speaker_scaled_time):VowelGOOSE 1.002 1.005 0.572 0.449558
## s(speaker_scaled_time):VowelKIT 2.260 2.817 2.301 0.104536
## s(speaker_scaled_time):VowelLOT 2.151 2.685 10.693 5.63e-06 ***
## s(speaker_scaled_time):VowelNURSE 1.003 1.005 0.874 0.348830
## s(speaker_scaled_time):VowelSTART 2.123 2.652 1.005 0.291813
## s(speaker_scaled_time):VowelSTRUT 5.784 6.938 3.088 0.002243 **
## s(speaker_scaled_time):VowelTHOUGHT 1.001 1.002 10.225 0.001377 **
## s(speaker_scaled_time):VowelTRAP 2.472 3.084 7.413 5.51e-05 ***
## s(speaker_scaled_art_rate):VowelDRESS 1.001 1.003 0.024 0.881028
## s(speaker_scaled_art_rate):VowelFLEECE 1.873 2.375 2.373 0.107352
## s(speaker_scaled_art_rate):VowelFOOT 1.774 2.239 6.239 0.000900 ***
## s(speaker_scaled_art_rate):VowelGOOSE 1.009 1.017 31.247 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelKIT 1.003 1.006 47.259 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelLOT 1.003 1.005 37.963 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelNURSE 1.001 1.003 0.054 0.818564
## s(speaker_scaled_art_rate):VowelSTART 1.003 1.006 28.939 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelSTRUT 1.003 1.006 33.869 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelTHOUGHT 1.003 1.006 10.960 0.000917 ***
## s(speaker_scaled_art_rate):VowelTRAP 1.025 1.050 15.676 5.54e-05 ***
## s(speaker_scaled_amp_max):VowelDRESS 1.013 1.025 43.623 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelFLEECE 2.730 3.469 34.134 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelFOOT 1.005 1.009 143.998 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelGOOSE 1.600 2.006 82.689 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelKIT 4.334 5.355 73.701 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelLOT 4.140 5.133 263.719 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelNURSE 1.003 1.006 41.461 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelSTART 3.346 4.201 26.595 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelSTRUT 3.983 4.945 248.067 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelTHOUGHT 1.801 2.290 23.570 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelTRAP 3.855 4.811 48.691 < 2e-16 ***
## s(speaker_scaled_pitch):VowelDRESS 4.562 5.602 9.063 < 2e-16 ***
## s(speaker_scaled_pitch):VowelFLEECE 3.448 4.327 5.337 0.000170 ***
## s(speaker_scaled_pitch):VowelFOOT 3.672 4.558 11.608 < 2e-16 ***
## s(speaker_scaled_pitch):VowelGOOSE 3.452 4.321 5.568 0.000162 ***
## s(speaker_scaled_pitch):VowelKIT 5.503 6.606 31.202 < 2e-16 ***
## s(speaker_scaled_pitch):VowelLOT 5.146 6.237 22.858 < 2e-16 ***
## s(speaker_scaled_pitch):VowelNURSE 4.068 5.038 8.869 < 2e-16 ***
## s(speaker_scaled_pitch):VowelSTART 1.004 1.009 2.657 0.101772
## s(speaker_scaled_pitch):VowelSTRUT 6.966 7.970 19.555 < 2e-16 ***
## s(speaker_scaled_pitch):VowelTHOUGHT 3.691 4.606 9.758 < 2e-16 ***
## s(speaker_scaled_pitch):VowelTRAP 4.180 5.169 5.877 1.98e-05 ***
## s(Speaker) 191.609 214.000 5065.674 < 2e-16 ***
## s(Vowel,Speaker) 1988.359 2361.000 82.424 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.755 Deviance explained = 61.9%
## fREML = 1.857e+05 Scale est. = 1 n = 109343
There are no obvious systematic differences in the model output. It explains somewhat less deviance, but this is not a cause for concern if the model still carries the main effects we are interested in and fits the data better.
With the t-distribution, quite a bit more is required in terms of degrees of
freedom to capture the vowel space of each speaker (see s(Speaker) and
s(Vowel, Speaker)).
The same items seem to appear as significant in our list of smooth terms as in Structure 1. If anything, there is a tendency for this model to indicate a higher degree of significance to these terms.
We will look at amplitude first, as it is our main variable of interest.
gamm_fit_preds <- get_predictions(
gamm_fit_scat,
cond = list(
'speaker_scaled_amp_max' = seq(-2.5, 2.5, length.out=200),
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
gamm_fit_preds %>%
mutate(# Enable plot to distinguish between high and low vowels (for faceting)
height = if_else(Vowel %in% high_vowels, "high", "low or mid")
) %>%
ggplot(
aes(
y = fit,
x = speaker_scaled_amp_max,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.3,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = "F1 and Amplitude by Vowel",
x = "Scaled amplitude",
y = "F1"
)
Figure 2.21: F1 and amplitude by vowel, scaled T model.
# facet_wrap(
# facets=vars(height), scales = "free"
# )
Figure 2.21 Is similar to the previous plots of the same sort. Again, lot and trap, seem to merge rather than cross. We will keep this in mind in interpreting the effect of amplitude on the vowel space as a whole.
We now run diagnostic checks on the model:Figure 2.22: GAM model checks.
The text output again shows that our choice of k is fine.
The QQ plot and histogram of residuals suggest that the performance of our residuals is superior.
The above models are fit with fREML. The resulting p-values are not
necessarily reliable. We use ML and the model comparison based significance
test recommended by Sóskuthy (2017).
We refit models using ML to perform significance testing by model comparison. This takes a long time to fit (around 40-50 minutes for each model).
We use structure 1, the model structure with random intercepts by speaker and by
speaker, by vowel. The random intercept structure can be thought of as
calibrating the model for the overall vowel space of each speaker. While
the previous section showed good evidence that a scaled t distribution would
be superior for our residuals, attempts to fit scaled t models using ML were
unsuccessful due to limitations of both time and memory.
gamm_fit_full <- bam(
F1_50 ~
participant_gender +
Vowel +
s(speaker_scaled_time, by=Vowel) +
s(speaker_scaled_art_rate, by=Vowel) +
s(speaker_scaled_amp_max, by=Vowel) +
s(speaker_scaled_pitch, by=Vowel) +
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="ML"
)
write_rds(gamm_fit_full, here('models', 'gamm_fit_full.rds'))
# We will look at the summary for the full model to check that it doesn't
# radically differ from the fREML model Structure 1.
gamm_fit_full_summary <- summary(gamm_fit_full)
write_rds(gamm_fit_full_summary, here('models', 'gamm_fit_full_summary.rds'))
gamm_fit_no_time <- bam(
F1_50 ~
participant_gender +
Vowel +
s(speaker_scaled_art_rate, by=Vowel) +
s(speaker_scaled_amp_max, by=Vowel) +
s(speaker_scaled_pitch, by=Vowel) +
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="ML"
)
write_rds(gamm_fit_no_time, here('models', 'gamm_fit_no_time.rds'))
gamm_fit_no_art <- bam(
F1_50 ~
participant_gender +
Vowel +
s(speaker_scaled_time, by=Vowel) +
s(speaker_scaled_amp_max, by=Vowel) +
s(speaker_scaled_pitch, by=Vowel) +
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="ML"
)
write_rds(gamm_fit_no_art, here('models', 'gamm_fit_no_art.rds'))
gamm_fit_no_amp <- bam(
F1_50 ~
participant_gender +
Vowel +
s(speaker_scaled_time, by=Vowel) +
s(speaker_scaled_art_rate, by=Vowel) +
s(speaker_scaled_pitch, by=Vowel) +
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="ML"
)
write_rds(gamm_fit_no_amp, here('models', 'gamm_fit_no_amp.rds'))
gamm_fit_no_pitch <- bam(
F1_50 ~
participant_gender +
Vowel +
s(speaker_scaled_time, by=Vowel) +
s(speaker_scaled_art_rate, by=Vowel) +
s(speaker_scaled_amp_max, by=Vowel) +
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="ML"
)
write_rds(gamm_fit_no_pitch, here('models', 'gamm_fit_no_pitch.rds'))
We generate p-values for each of the variables using model comparison.
gamm_fit_full <- read_rds(here('models', 'gamm_fit_full.rds'))
gamm_fit_full_summary <- read_rds(here('models', 'gamm_fit_full_summary.rds'))
# time.
gamm_fit_no_time <- read_rds(here('models', 'gamm_fit_no_time.rds'))
compareML(gamm_fit_full, gamm_fit_no_time)
## gamm_fit_full: F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(Speaker,
## bs = "re") + s(Speaker, Vowel, bs = "re")
##
## gamm_fit_no_time: F1_50 ~ participant_gender + Vowel + s(speaker_scaled_art_rate,
## by = Vowel) + s(speaker_scaled_amp_max, by = Vowel) + s(speaker_scaled_pitch,
## by = Vowel) + s(Speaker, bs = "re") + s(Speaker, Vowel, bs = "re")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 gamm_fit_no_time 624671.0 80
## 2 gamm_fit_full 624634.5 102 36.498 22.000 2.215e-07 ***
##
## AIC difference: -58.29, model gamm_fit_full has lower AIC.
# articulation rate.
gamm_fit_no_art <- read_rds(here('models', 'gamm_fit_no_art.rds'))
compareML(gamm_fit_full, gamm_fit_no_art)
## gamm_fit_full: F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(Speaker,
## bs = "re") + s(Speaker, Vowel, bs = "re")
##
## gamm_fit_no_art: F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_amp_max, by = Vowel) + s(speaker_scaled_pitch,
## by = Vowel) + s(Speaker, bs = "re") + s(Speaker, Vowel, bs = "re")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 gamm_fit_no_art 624713.1 80
## 2 gamm_fit_full 624634.5 102 78.602 22.000 < 2e-16 ***
##
## AIC difference: -132.77, model gamm_fit_full has lower AIC.
# amplitude.
gamm_fit_no_amp <- read_rds(here('models', 'gamm_fit_no_amp.rds'))
compareML(gamm_fit_full, gamm_fit_no_amp)
## gamm_fit_full: F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(Speaker,
## bs = "re") + s(Speaker, Vowel, bs = "re")
##
## gamm_fit_no_amp: F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_pitch,
## by = Vowel) + s(Speaker, bs = "re") + s(Speaker, Vowel, bs = "re")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 gamm_fit_no_amp 626422.7 80
## 2 gamm_fit_full 624634.5 102 1788.190 22.000 < 2e-16 ***
##
## AIC difference: -3614.18, model gamm_fit_full has lower AIC.
# pitch
gamm_fit_no_pitch <- read_rds(here('models', 'gamm_fit_no_pitch.rds'))
compareML(gamm_fit_full, gamm_fit_no_pitch)
## gamm_fit_full: F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(Speaker,
## bs = "re") + s(Speaker, Vowel, bs = "re")
##
## gamm_fit_no_pitch: F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(Speaker, bs = "re") + s(Speaker, Vowel, bs = "re")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 gamm_fit_no_pitch 624893.4 80
## 2 gamm_fit_full 624634.5 102 258.941 22.000 < 2e-16 ***
##
## AIC difference: -566.26, model gamm_fit_full has lower AIC.
The four tables above present the contribution of each of our variables to the model. Adding new variables will always result in a better fit to the data. The question is whether the increase in model complexity is worth it for the increase in information about our variable of pseudointerest. For, for instance, our time variable requires, in effect, around 22 additional parameters over the model without time. The information given here is taken by the \(\Chi^2\) test to be significant.
The same is true for our other variable. Note that the Difference column can
be used to estimate the respective magnitude of the contribution of each
variable, with amplitude providing the most information for the cost in degrees
of freedom, and time providing the least.
We use the model fit with ML to generate plots of our model predictions for each vowel and for each of our main fixed variables.
We first define a function to plot predictions from the models and to thereby reduce repetition of code.
all_vowel_plot <- function(predictions, xvar, xlabel, plot_title) {
xvar <- enquo(xvar)
predictions %>%
mutate(# Enable plot to distinguish between high and low vowels (for faceting)
height = if_else(Vowel %in% high_vowels, "high", "low or mid")
) %>%
ggplot(
aes(
y = fit,
x = !!xvar,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.3,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = plot_title,
x = xlabel,
y = "F1"
)
# facet_wrap(
# facets=vars(height), scales = "free"
# )
}
The following figure is Figure 3 in the paper.
gamm_fit_preds <- get_predictions(
gamm_fit_full,
cond = list(
'speaker_scaled_amp_max' = seq(-2.5, 2.5, length.out=200),
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Speaker,Vowel)
##
figure_3 <- all_vowel_plot(
gamm_fit_preds,
speaker_scaled_amp_max,
"Scaled amplitude",
"F1 and Amplitude by Vowel"
)
ggsave(
here('plots', 'figure_3.png'),
plot = figure_3,
dpi = 500,
units = "cm",
width = 20,
height = 12.5
)
figure_3
Figure 2.23: F1 and amplitude by vowel, ML model.
Figure @red(fig:gamm-full-amp-plot) shows strong evidence for increases in F1 accompanying increases in amplitude. This effect is stronger for some vowels than for others. The estimates for lot and trap cross one another as amplitude increases. We also see that some clusters of vowels are distinguishable in F1 at some amplitudes and not at others. These effects are concerning because they would lead an analyst to reach different conclusions about the configuration of the vowel space at higher versus lower amplitudes.
gamm_fit_preds <- get_predictions(
gamm_fit_full,
cond = list(
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = seq(-2.5, 2.5, length.out=200),
'speaker_scaled_time' = 0.5,
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Speaker,Vowel)
##
all_vowel_plot(gamm_fit_preds, speaker_scaled_art_rate, "Scaled articulation rate", "F1 and Articulation Rate by Vowel")
Figure 2.24: F1 and articulation rate by vowel.
In Figure 2.24 the main thing we are capturing seems to be a reduction in vowel space, with thought and goose moving in the opposite direction to the other vowels. That is, we have contraction occurring at higher rates of speech. Note that the magnitude of this effect is much smaller than that of amplitude (Figure @ref(fig:gamm_full_amp_plot)).
gamm_fit_preds <- get_predictions(
gamm_fit_full,
cond = list(
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F',
'speaker_scaled_pitch' = seq(-2.5, 2.5, length.out=200),
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Speaker,Vowel)
##
all_vowel_plot(gamm_fit_preds, speaker_scaled_pitch, "Scaled pitch", "F1 and Pitch by Vowel")
Figure 2.25: F1 and pitch by vowel.
There is an interesting U shaped curve for many vowels when we look at pitch.
gamm_fit_preds <- get_predictions(
gamm_fit_full,
cond = list(
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = seq(0, 1, length.out=200),
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Speaker,Vowel)
##
all_vowel_plot(gamm_fit_preds, speaker_scaled_time, "Scaled time", "F1 and Time by Vowel")
Figure 2.26: F1 and time by vowel.
It is very hard to take this plot as presenting any good evidence for an effect of time through a monologue and F1. Certainly, it is the variable our models are least confident about. If anything, there is a small effect where lot and trap reducing over the course of the monologue and kit increasing.
We now look at our GAMM diagnostics.
gam.check(gamm_fit_full)
##
## Method: ML Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-0.3467237,0.1643703]
## (score 624634.5 & scale 5063.495).
## Hessian positive definite, eigenvalue range [0.0174575,54689.18].
## Model rank = 3000 / 3000
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time):VowelDRESS 9.00 1.02 1.00 0.56
## s(speaker_scaled_time):VowelFLEECE 9.00 1.02 1.00 0.56
## s(speaker_scaled_time):VowelFOOT 9.00 1.03 1.00 0.66
## s(speaker_scaled_time):VowelGOOSE 9.00 1.03 1.00 0.56
## s(speaker_scaled_time):VowelKIT 9.00 1.04 1.00 0.58
## s(speaker_scaled_time):VowelLOT 9.00 1.81 1.00 0.59
## s(speaker_scaled_time):VowelNURSE 9.00 1.03 1.00 0.59
## s(speaker_scaled_time):VowelSTART 9.00 1.03 1.00 0.60
## s(speaker_scaled_time):VowelSTRUT 9.00 5.93 1.00 0.59
## s(speaker_scaled_time):VowelTHOUGHT 9.00 1.02 1.00 0.58
## s(speaker_scaled_time):VowelTRAP 9.00 1.96 1.00 0.58
## s(speaker_scaled_art_rate):VowelDRESS 9.00 1.02 1.00 0.55
## s(speaker_scaled_art_rate):VowelFLEECE 9.00 1.03 1.00 0.61
## s(speaker_scaled_art_rate):VowelFOOT 9.00 1.03 1.00 0.56
## s(speaker_scaled_art_rate):VowelGOOSE 9.00 1.02 1.00 0.55
## s(speaker_scaled_art_rate):VowelKIT 9.00 1.03 1.00 0.56
## s(speaker_scaled_art_rate):VowelLOT 9.00 1.02 1.00 0.59
## s(speaker_scaled_art_rate):VowelNURSE 9.00 1.02 1.00 0.54
## s(speaker_scaled_art_rate):VowelSTART 9.00 1.02 1.00 0.59
## s(speaker_scaled_art_rate):VowelSTRUT 9.00 2.19 1.00 0.57
## s(speaker_scaled_art_rate):VowelTHOUGHT 9.00 1.02 1.00 0.56
## s(speaker_scaled_art_rate):VowelTRAP 9.00 1.03 1.00 0.60
## s(speaker_scaled_amp_max):VowelDRESS 9.00 1.03 1.00 0.58
## s(speaker_scaled_amp_max):VowelFLEECE 9.00 1.05 1.00 0.56
## s(speaker_scaled_amp_max):VowelFOOT 9.00 1.02 1.00 0.58
## s(speaker_scaled_amp_max):VowelGOOSE 9.00 1.04 1.00 0.54
## s(speaker_scaled_amp_max):VowelKIT 9.00 3.50 1.00 0.47
## s(speaker_scaled_amp_max):VowelLOT 9.00 3.74 1.00 0.58
## s(speaker_scaled_amp_max):VowelNURSE 9.00 1.02 1.00 0.57
## s(speaker_scaled_amp_max):VowelSTART 9.00 1.04 1.00 0.51
## s(speaker_scaled_amp_max):VowelSTRUT 9.00 5.06 1.00 0.56
## s(speaker_scaled_amp_max):VowelTHOUGHT 9.00 1.02 1.00 0.55
## s(speaker_scaled_amp_max):VowelTRAP 9.00 2.63 1.00 0.54
## s(speaker_scaled_pitch):VowelDRESS 9.00 3.53 1.03 0.98
## s(speaker_scaled_pitch):VowelFLEECE 9.00 1.02 1.03 0.97
## s(speaker_scaled_pitch):VowelFOOT 9.00 3.00 1.03 0.97
## s(speaker_scaled_pitch):VowelGOOSE 9.00 2.70 1.03 0.97
## s(speaker_scaled_pitch):VowelKIT 9.00 6.06 1.03 0.98
## s(speaker_scaled_pitch):VowelLOT 9.00 4.67 1.03 0.98
## s(speaker_scaled_pitch):VowelNURSE 9.00 2.86 1.03 0.95
## s(speaker_scaled_pitch):VowelSTART 9.00 1.02 1.03 1.00
## s(speaker_scaled_pitch):VowelSTRUT 9.00 7.55 1.03 1.00
## s(speaker_scaled_pitch):VowelTHOUGHT 9.00 2.81 1.03 0.97
## s(speaker_scaled_pitch):VowelTRAP 9.00 3.02 1.03 0.98
## s(Speaker) 216.00 192.06 NA NA
## s(Speaker,Vowel) 2376.00 1912.84 NA NA
All basically equivalent to our fREML model (Structure 1) as presented above.
gamm_fit_full_summary
##
## Family: gaussian
## Link function: identity
##
## Formula:
## F1_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(Speaker,
## bs = "re") + s(Speaker, Vowel, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 432.283 4.036 107.096 < 2e-16 ***
## participant_genderM -55.258 5.608 -9.854 < 2e-16 ***
## VowelFLEECE -10.894 3.785 -2.878 0.00400 **
## VowelFOOT 59.716 4.051 14.741 < 2e-16 ***
## VowelGOOSE -22.286 3.872 -5.755 8.67e-09 ***
## VowelKIT 82.505 3.789 21.774 < 2e-16 ***
## VowelLOT 159.218 3.819 41.686 < 2e-16 ***
## VowelNURSE 11.596 3.912 2.964 0.00304 **
## VowelSTART 358.035 3.892 91.989 < 2e-16 ***
## VowelSTRUT 255.220 3.790 67.339 < 2e-16 ***
## VowelTHOUGHT 9.275 3.831 2.421 0.01547 *
## VowelTRAP 145.829 3.811 38.270 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(speaker_scaled_time):VowelDRESS 1.023 1.046 6.288 0.010285 *
## s(speaker_scaled_time):VowelFLEECE 1.023 1.046 2.519 0.108771
## s(speaker_scaled_time):VowelFOOT 1.030 1.060 0.905 0.334325
## s(speaker_scaled_time):VowelGOOSE 1.026 1.051 0.000 0.999916
## s(speaker_scaled_time):VowelKIT 1.043 1.084 3.769 0.044233 *
## s(speaker_scaled_time):VowelLOT 1.813 2.264 12.395 3.70e-06 ***
## s(speaker_scaled_time):VowelNURSE 1.031 1.060 0.533 0.479569
## s(speaker_scaled_time):VowelSTART 1.029 1.057 0.014 0.953339
## s(speaker_scaled_time):VowelSTRUT 5.929 7.087 4.270 0.000145 ***
## s(speaker_scaled_time):VowelTHOUGHT 1.022 1.043 5.807 0.014659 *
## s(speaker_scaled_time):VowelTRAP 1.959 2.448 8.387 0.000157 ***
## s(speaker_scaled_art_rate):VowelDRESS 1.025 1.049 0.021 0.942557
## s(speaker_scaled_art_rate):VowelFLEECE 1.034 1.067 2.255 0.138682
## s(speaker_scaled_art_rate):VowelFOOT 1.026 1.051 6.749 0.008898 **
## s(speaker_scaled_art_rate):VowelGOOSE 1.024 1.048 22.923 1.14e-06 ***
## s(speaker_scaled_art_rate):VowelKIT 1.027 1.053 37.875 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelLOT 1.016 1.031 25.116 3.94e-07 ***
## s(speaker_scaled_art_rate):VowelNURSE 1.024 1.048 0.155 0.721426
## s(speaker_scaled_art_rate):VowelSTART 1.020 1.040 18.161 1.70e-05 ***
## s(speaker_scaled_art_rate):VowelSTRUT 2.189 2.777 9.581 1.00e-05 ***
## s(speaker_scaled_art_rate):VowelTHOUGHT 1.023 1.046 5.623 0.016229 *
## s(speaker_scaled_art_rate):VowelTRAP 1.030 1.059 10.693 0.000992 ***
## s(speaker_scaled_amp_max):VowelDRESS 1.030 1.060 23.348 9.02e-07 ***
## s(speaker_scaled_amp_max):VowelFLEECE 1.054 1.106 53.689 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelFOOT 1.017 1.033 96.898 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelGOOSE 1.037 1.073 91.535 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelKIT 3.496 4.398 57.837 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelLOT 3.738 4.686 279.917 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelNURSE 1.024 1.048 21.159 3.28e-06 ***
## s(speaker_scaled_amp_max):VowelSTART 1.037 1.072 87.429 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelSTRUT 5.057 6.206 238.074 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelTHOUGHT 1.018 1.037 28.550 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelTRAP 2.625 3.340 62.677 < 2e-16 ***
## s(speaker_scaled_pitch):VowelDRESS 3.531 4.436 5.675 0.000113 ***
## s(speaker_scaled_pitch):VowelFLEECE 1.020 1.041 9.637 0.001884 **
## s(speaker_scaled_pitch):VowelFOOT 3.000 3.773 8.876 5.46e-07 ***
## s(speaker_scaled_pitch):VowelGOOSE 2.701 3.423 3.657 0.009247 **
## s(speaker_scaled_pitch):VowelKIT 6.060 7.247 24.955 < 2e-16 ***
## s(speaker_scaled_pitch):VowelLOT 4.671 5.776 23.422 < 2e-16 ***
## s(speaker_scaled_pitch):VowelNURSE 2.862 3.609 5.534 0.000522 ***
## s(speaker_scaled_pitch):VowelSTART 1.024 1.047 3.000 0.080781 .
## s(speaker_scaled_pitch):VowelSTRUT 7.545 8.475 19.149 < 2e-16 ***
## s(speaker_scaled_pitch):VowelTHOUGHT 2.808 3.553 5.253 0.000739 ***
## s(speaker_scaled_pitch):VowelTRAP 3.022 3.822 3.684 0.006279 **
## s(Speaker) 192.059 214.000 2231.314 < 2e-16 ***
## s(Speaker,Vowel) 1912.840 2361.000 37.603 4.89e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.759 Deviance explained = 76.4%
## -ML = 6.2463e+05 Scale est. = 5063.5 n = 109343
The summary values are also roughly equivalent.
We now fit GAMM models to the data from each vowel independently to see if this has any effect of our conclusions. For these modules, it is computationally feasible to fit with random smooths for each speaker.
Again, we set this block to not evaluate and load precalculated results in the following cell.
vowel_gam_models <- qb_vowels %>%
# Group by vowel and nest to create a column of dataframes corresponding
# to each vowel.
group_by(Vowel) %>%
nest() %>%
# Apply gam model.
mutate(
model = map(
data,
~ bam(
F1_50 ~
participant_gender +
s(speaker_scaled_time) +
s(speaker_scaled_art_rate) +
s(speaker_scaled_amp_max) +
s(speaker_scaled_pitch) +
s(speaker_scaled_time, Speaker, bs="fs", k=5, m=1) + #k=5 to reduce load.
s(Speaker, bs="re"),
data = .x,
method="fREML",
discrete = TRUE,
nthreads = 8 # Change this depending on number of cores available.
)
)
)
write_rds(vowel_gam_models, here('models', 'vowel_gam_models.rds'))
vowel_gam_models <- read_rds(here('models', 'vowel_gam_models.rds'))
We first look to see if the residual distribution looks any better.
walk2(
vowel_gam_models$Vowel,
vowel_gam_models$model,
~ qq.gam(.y, main=glue('{.x}: Resids vs. linear pred.'))
)
It is not unexpected that each of these vowels has a distribution which has
strange behaviour in the tails. Some vowels do not have much room to move ‘up’
and some do not have much room to move ‘down’. So, for instance, there is
more variation possible at low F1 values for
start and
strut, and the model
thus seems to perform more poorly at low values.
If we look at the model for start in more detail we see:
gam.check(
vowel_gam_models %>%
filter(
Vowel == "START"
) %>%
pull(model) %>%
pluck(1)
)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] -6.173218e-05 -4.544179e-05 -9.157447e-05 -4.036669e-05 -1.841179e-07
## [6] -5.416620e-06 -5.416617e-06 1.249726e-04
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## 6.172205e-05 4.976852e-12 -2.294380e-10 -3.782817e-11 -4.798869e-06
## 4.976852e-12 4.543877e-05 8.685806e-12 -4.024387e-12 -1.033302e-06
## -2.294380e-10 8.685806e-12 9.156081e-05 5.810362e-11 5.979093e-06
## -3.782817e-11 -4.024387e-12 5.810362e-11 4.036320e-05 -9.930338e-07
## -4.798869e-06 -1.033302e-06 5.979093e-06 -9.930338e-07 8.079979e+00
## -1.662477e-06 -5.245999e-07 8.928690e-06 -1.845007e-06 -1.000281e-01
## -1.662477e-06 -5.245999e-07 8.928690e-06 -1.845007e-06 -1.000281e-01
## d -3.772633e-05 -3.269180e-06 -7.169194e-05 -1.549775e-05 -3.931401e+01
## [,6] [,7] [,8]
## -1.662477e-06 -1.662477e-06 -3.772633e-05
## -5.245999e-07 -5.245999e-07 -3.269180e-06
## 8.928690e-06 8.928690e-06 -7.169194e-05
## -1.845007e-06 -1.845007e-06 -1.549775e-05
## -1.000281e-01 -1.000281e-01 -3.931401e+01
## 2.281001e+01 2.281001e+01 -4.974211e+01
## 2.281001e+01 2.281001e+01 -4.974211e+01
## d -4.974211e+01 -4.974211e+01 2.933500e+03
##
## Model rank = 1334 / 1334
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.0 1.0 1.02 0.93
## s(speaker_scaled_art_rate) 9.0 1.0 0.99 0.20
## s(speaker_scaled_amp_max) 9.0 1.0 1.01 0.81
## s(speaker_scaled_pitch) 9.0 1.0 0.99 0.28
## s(speaker_scaled_time,Speaker) 1080.0 178.1 1.02 0.89
## s(Speaker) 216.0 99.5 NA NA
There does seem to be something strange happening for some speakers, where the model predicts low F1 values, but the actual value is even lower. Again, the general pattern looks OK. We are talking about a handful of affected observations.
We check our k values.
vowel_gam_models <- vowel_gam_models %>%
mutate(
kcheck = map(model, k.check)
)
print_kcheck <- function(vowel, kcheck) { # Taken from mgcv gam.check code.
cat(glue('{vowel}\n'))
cat("\nBasis dimension (k) checking results. Low p-value (k-index<1) may\n")
cat("indicate that k is too low, especially if edf is close to k\'.\n\n")
printCoefmat(kcheck, digits=3)
cat('\n')
}
walk2(
as.character(vowel_gam_models$Vowel),
vowel_gam_models$kcheck,
print_kcheck
)
## DRESS
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.0 1.9 1.01 0.74
## s(speaker_scaled_art_rate) 9.0 1.0 0.99 0.35
## s(speaker_scaled_amp_max) 9.0 1.0 1.00 0.48
## s(speaker_scaled_pitch) 9.0 5.0 0.98 0.13
## s(speaker_scaled_time,Speaker) 1080.0 321.9 1.01 0.73
## s(Speaker) 216.0 104.4 NA NA
##
## NURSE
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 1.00 1.00 0.54
## s(speaker_scaled_art_rate) 9.00 1.00 1.01 0.83
## s(speaker_scaled_amp_max) 9.00 1.00 1.00 0.38
## s(speaker_scaled_pitch) 9.00 5.06 0.98 0.05 *
## s(speaker_scaled_time,Speaker) 1080.00 200.52 1.00 0.58
## s(Speaker) 216.00 103.07 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## KIT
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 4.76 0.99 0.2600
## s(speaker_scaled_art_rate) 9.00 1.00 1.00 0.4775
## s(speaker_scaled_amp_max) 9.00 3.54 0.98 0.0975 .
## s(speaker_scaled_pitch) 9.00 5.50 0.96 0.0075 **
## s(speaker_scaled_time,Speaker) 1080.00 237.57 0.99 0.2650
## s(Speaker) 216.00 99.88 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## TRAP
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 2.32 0.99 0.365
## s(speaker_scaled_art_rate) 9.00 1.00 0.99 0.185
## s(speaker_scaled_amp_max) 9.00 3.80 1.01 0.743
## s(speaker_scaled_pitch) 9.00 4.71 0.98 0.058 .
## s(speaker_scaled_time,Speaker) 1080.00 263.14 0.99 0.378
## s(Speaker) 216.00 102.56 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## LOT
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 2.27 0.99 0.13
## s(speaker_scaled_art_rate) 9.00 1.00 0.99 0.36
## s(speaker_scaled_amp_max) 9.00 3.83 0.99 0.34
## s(speaker_scaled_pitch) 9.00 4.75 1.03 0.98
## s(speaker_scaled_time,Speaker) 1080.00 184.03 0.99 0.14
## s(Speaker) 216.00 99.27 NA NA
##
## THOUGHT
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 1.00 0.98 0.060 .
## s(speaker_scaled_art_rate) 9.00 1.00 0.98 0.138
## s(speaker_scaled_amp_max) 9.00 2.48 1.00 0.610
## s(speaker_scaled_pitch) 9.00 3.72 0.99 0.177
## s(speaker_scaled_time,Speaker) 1080.00 189.96 0.98 0.043 *
## s(Speaker) 216.00 100.01 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## FLEECE
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 1.66 0.98 0.068 .
## s(speaker_scaled_art_rate) 9.00 1.41 0.99 0.300
## s(speaker_scaled_amp_max) 9.00 3.03 1.02 0.922
## s(speaker_scaled_pitch) 9.00 3.81 0.99 0.210
## s(speaker_scaled_time,Speaker) 1080.00 251.93 0.98 0.098 .
## s(Speaker) 216.00 102.28 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## STRUT
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 1.83 0.99 0.28
## s(speaker_scaled_art_rate) 9.00 1.65 1.03 0.93
## s(speaker_scaled_amp_max) 9.00 3.20 0.98 0.14
## s(speaker_scaled_pitch) 9.00 5.52 0.96 0.01 **
## s(speaker_scaled_time,Speaker) 1080.00 233.68 0.99 0.28
## s(Speaker) 216.00 100.57 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## START
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.0 1.0 1.03 0.96
## s(speaker_scaled_art_rate) 9.0 1.0 0.98 0.06 .
## s(speaker_scaled_amp_max) 9.0 1.0 1.01 0.80
## s(speaker_scaled_pitch) 9.0 1.0 1.00 0.57
## s(speaker_scaled_time,Speaker) 1080.0 178.1 1.03 0.95
## s(Speaker) 216.0 99.5 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## GOOSE
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 1.00 1.00 0.365
## s(speaker_scaled_art_rate) 9.00 2.29 0.98 0.077 .
## s(speaker_scaled_amp_max) 9.00 3.68 1.00 0.485
## s(speaker_scaled_pitch) 9.00 3.48 0.99 0.193
## s(speaker_scaled_time,Speaker) 1080.00 203.60 1.00 0.403
## s(Speaker) 216.00 92.97 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## FOOT
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(speaker_scaled_time) 9.00 1.00 1.00 0.45
## s(speaker_scaled_art_rate) 9.00 1.00 0.99 0.34
## s(speaker_scaled_amp_max) 9.00 3.26 0.99 0.28
## s(speaker_scaled_pitch) 9.00 3.58 1.01 0.78
## s(speaker_scaled_time,Speaker) 1065.00 136.00 1.00 0.48
## s(Speaker) 213.00 84.13 NA NA
All look OK!
We now plot the results together.
all_vowel_plots <- vowel_gam_models %>%
mutate(
preds = map(
model,
~ get_predictions(
.x,
cond = list(
'speaker_scaled_amp_max' = seq(-2.5, 2.5, length.out=200),
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'participant_gender' = 'F'
)
)
)
) %>%
select(Vowel, preds) %>%
unnest() %>%
mutate(
height = if_else(Vowel %in% high_vowels, 'high', 'mid or low')
) %>%
ggplot(
aes(
y = fit,
x = speaker_scaled_amp_max,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.25,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = "F1 and Amplitude by Vowel",
x = "Scaled amplitude",
y = "F1"
) +
facet_wrap(
facets=vars(height), scales = "free"
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_F_408. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_384. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_753. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_F_408. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
all_vowel_plots
Figure 2.27: Amplitude predictions from by-vowel models.
There is a more extreme reduction in F1 values for goose at low values of scaled amplitude. The values of lot and trap get a little closer together at low values as well. However, the overall story is the same. We see increase in F1 with amplitude.
all_vowel_plots <- vowel_gam_models %>%
mutate(
preds = map(
model,
~ get_predictions(
.x,
cond = list(
'speaker_scaled_pitch' = seq(-2.5, 2.5, length.out=200),
'speaker_scaled_amp_max' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'participant_gender' = 'F'
)
)
)
) %>%
select(Vowel, preds) %>%
unnest() %>%
mutate(
height = if_else(Vowel %in% high_vowels, 'high', 'mid or low')
) %>%
ggplot(
aes(
y = fit,
x = speaker_scaled_pitch,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.25,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = "F1 and Pitch by Vowel",
x = "Scaled pitch",
y = "F1"
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_F_408. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_384. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_753. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_F_408. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
# facet_wrap(
# facets=vars(height), scales = "free"
# )
all_vowel_plots
Figure 2.28: Pitch predictions from by-vowel models.
The U shaped plots again.
all_vowel_plots <- vowel_gam_models %>%
mutate(
preds = map(
model,
~ get_predictions(
.x,
cond = list(
'speaker_scaled_art_rate' = seq(-2.5, 2.5, length.out=200),
'speaker_scaled_amp_max' = 0,
'speaker_scaled_pitch' = 0,
'speaker_scaled_time' = 0.5,
'participant_gender' = 'F'
)
)
)
) %>%
select(Vowel, preds) %>%
unnest() %>%
mutate(
height = if_else(Vowel %in% high_vowels, 'high', 'mid or low')
) %>%
ggplot(
aes(
y = fit,
x = speaker_scaled_art_rate,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.25,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = "F1 and Articulation Rate by Vowel",
x = "Scaled articulation rate",
y = "F1"
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_F_408. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_384. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_753. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_F_408. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
# facet_wrap(
# facets=vars(height), scales = "free"
# )
all_vowel_plots
Figure 2.29: Articulation rate predictions from by-vowel models.
Interestingly, if we model the vowels independently, we lose the increase in F1 with articulation rate of thought. However, the rest of the phenoemena explicable by vowel space contraction and expansion with articulation rate remaine (Figure 2.29).
all_vowel_plots <- vowel_gam_models %>%
mutate(
preds = map(
model,
~ get_predictions(
.x,
cond = list(
'speaker_scaled_time' = seq(0, 1, length.out=200),
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F'
)
)
)
) %>%
select(Vowel, preds) %>%
unnest() %>%
mutate(
height = if_else(Vowel %in% high_vowels, 'high', 'mid or low')
) %>%
ggplot(
aes(
y = fit,
x = speaker_scaled_time,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.25,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = "F1 and Time by Vowel",
x = "Scaled time",
y = "F1"
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_F_408. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_384. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_753. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_F_408. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
all_vowel_plots
Figure 2.30: Time predictions from by-vowel models.
There is no great change in our predictions here, with the possible exception of an acceleration of change in kit, trap, and lot at the end of a monologue.
As we did with our simple linear models, we do some simple p-value significance testing to determine how often our various predictors come out as significant. This takes some time as it requires a summary to be generated for each gamm.
vowel_gam_sig_coeffs <- vowel_gam_models %>%
mutate(
summary = map(model, summary),
smooth_coefficients = map(
summary,
~ as_tibble(.x$s.table, rownames = 'variable'),
),
significant_variables = map(
smooth_coefficients,
~ .x %>% filter(`p-value` <= 0.05)
)
) %>%
select(
Vowel, significant_variables
) %>%
unnest(significant_variables)
write_rds(vowel_gam_sig_coeffs, here("models", "vowel_gam_sig_coeffs.rds"))
vowel_gam_sig_coeffs <- read_rds(here("models", "vowel_gam_sig_coeffs.rds"))
Let’s have a look at the values for speaker_scaled_amp_max.
vowel_gam_sig_coeffs %>%
filter(
variable == "s(speaker_scaled_amp_max)"
)
All 11 vowel have significant p-values for amplitude.
We compare these to the values for speaker_scaled_pitch:
vowel_gam_sig_coeffs %>%
filter(
variable == "s(speaker_scaled_pitch)"
)
All but start for pitch.
Now speaker_scaled_art_rate:
vowel_gam_sig_coeffs %>%
filter(
variable == "s(speaker_scaled_art_rate)"
)
All but dress, nurse, and fleece for articulation rate.
Finally, we look at s(speaker_scaled_time) and s(speaker_scaled_time, Speaker).
vowel_gam_sig_coeffs %>%
filter(
variable == 's(speaker_scaled_time)'
)
Only three vowels come out as significant for speaker_scaled_time.
Let’s have a look at these speaker_scaled_time smooths again:
all_vowel_plots <- vowel_gam_models %>%
filter(
Vowel %in% c('KIT', 'TRAP', 'STRUT', 'LOT')
) %>%
mutate(
preds = map(
model,
~ get_predictions(
.x,
cond = list(
'speaker_scaled_time' = seq(0, 1, length.out=200),
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F'
)
)
)
) %>%
select(Vowel, preds) %>%
unnest() %>%
ggplot(
aes(
y = fit,
x = speaker_scaled_time,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.25,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_with_foot
) +
scale_fill_manual(
values = vowel_colours_with_foot
) +
labs(
title = "F1 and Time by Vowel",
x = "Scaled time",
y = "F1"
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * speaker_scaled_time : numeric predictor; with 200 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speaker_scaled_time,Speaker),s(Speaker)
##
all_vowel_plots
Figure 2.31: Time predictions from by-vowel models.
All of these could have a line at a constant F1 value drawn through their confidence intervals. Whatever this effect is, it is very small compared to the effects of amplitude, pitch, and articulation rate.
We want to look at the consequences of amplitude change for the whole vowel space. In order to do this, we fit an fREML model for F2. For this purpose, we will match Structure 1.
gamm_fit_f2 <- bam(
F2_50 ~
participant_gender +
Vowel +
s(speaker_scaled_time, by=Vowel) +
s(speaker_scaled_art_rate, by=Vowel) +
s(speaker_scaled_amp_max, by=Vowel) +
s(speaker_scaled_pitch, by=Vowel) +
s(Speaker, bs="re") +
s(Speaker, Vowel, bs="re"),
data = qb_vowels,
method="fREML",
discrete = TRUE,
nthreads = 8 # Change this depending on number of cores available.
)
write_rds(gamm_fit_f2, here('models', 'gamm_fit_f2.rds'))
gamm_fit_f2_summary <- summary(gamm_fit_f2)
write_rds(gamm_fit_f2_summary, here('models', 'gamm_fit_f2_summary.rds'))
We look at the model summary.
gamm_fit_f2 <- read_rds(here('models', 'gamm_fit_f2.rds'))
gamm_fit_f2_summary <- read_rds(here('models', 'gamm_fit_f2_summary.rds'))
gamm_fit_f2_summary
##
## Family: gaussian
## Link function: identity
##
## Formula:
## F2_50 ~ participant_gender + Vowel + s(speaker_scaled_time, by = Vowel) +
## s(speaker_scaled_art_rate, by = Vowel) + s(speaker_scaled_amp_max,
## by = Vowel) + s(speaker_scaled_pitch, by = Vowel) + s(Speaker,
## bs = "re") + s(Speaker, Vowel, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2262.05 10.34 218.746 <2e-16 ***
## participant_genderM -238.98 12.00 -19.912 <2e-16 ***
## VowelFLEECE 14.81 11.70 1.266 0.205
## VowelFOOT -923.06 12.50 -73.832 <2e-16 ***
## VowelGOOSE -349.73 11.80 -29.628 <2e-16 ***
## VowelKIT -499.08 11.74 -42.496 <2e-16 ***
## VowelLOT -1059.53 11.98 -88.418 <2e-16 ***
## VowelNURSE -362.10 12.14 -29.820 <2e-16 ***
## VowelSTART -703.58 11.90 -59.109 <2e-16 ***
## VowelSTRUT -743.21 11.79 -63.061 <2e-16 ***
## VowelTHOUGHT -1254.25 12.41 -101.047 <2e-16 ***
## VowelTRAP -170.59 12.00 -14.218 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(speaker_scaled_time):VowelDRESS 5.096 6.205 3.687 0.001049 **
## s(speaker_scaled_time):VowelFLEECE 2.367 2.951 3.240 0.024344 *
## s(speaker_scaled_time):VowelFOOT 1.001 1.002 0.153 0.696292
## s(speaker_scaled_time):VowelGOOSE 2.979 3.703 3.388 0.010423 *
## s(speaker_scaled_time):VowelKIT 1.556 1.921 0.365 0.659576
## s(speaker_scaled_time):VowelLOT 1.001 1.002 1.642 0.199922
## s(speaker_scaled_time):VowelNURSE 1.441 1.754 1.289 0.191589
## s(speaker_scaled_time):VowelSTART 1.001 1.002 0.036 0.851360
## s(speaker_scaled_time):VowelSTRUT 1.001 1.001 3.589 0.058093 .
## s(speaker_scaled_time):VowelTHOUGHT 4.811 5.882 2.435 0.021576 *
## s(speaker_scaled_time):VowelTRAP 1.002 1.003 3.902 0.048212 *
## s(speaker_scaled_art_rate):VowelDRESS 1.959 2.483 26.449 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelFLEECE 1.005 1.010 30.132 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelFOOT 1.001 1.001 2.661 0.102717
## s(speaker_scaled_art_rate):VowelGOOSE 1.001 1.001 0.111 0.739731
## s(speaker_scaled_art_rate):VowelKIT 1.005 1.010 0.362 0.546487
## s(speaker_scaled_art_rate):VowelLOT 1.004 1.008 6.990 0.008004 **
## s(speaker_scaled_art_rate):VowelNURSE 1.967 2.487 2.311 0.093336 .
## s(speaker_scaled_art_rate):VowelSTART 1.001 1.002 3.988 0.045760 *
## s(speaker_scaled_art_rate):VowelSTRUT 1.001 1.002 5.908 0.015030 *
## s(speaker_scaled_art_rate):VowelTHOUGHT 4.328 5.346 8.648 < 2e-16 ***
## s(speaker_scaled_art_rate):VowelTRAP 1.023 1.046 0.010 0.977832
## s(speaker_scaled_amp_max):VowelDRESS 2.280 2.905 1.899 0.123316
## s(speaker_scaled_amp_max):VowelFLEECE 2.577 3.270 2.673 0.039636 *
## s(speaker_scaled_amp_max):VowelFOOT 1.004 1.008 0.355 0.555180
## s(speaker_scaled_amp_max):VowelGOOSE 1.001 1.003 132.986 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelKIT 1.001 1.003 1.222 0.269230
## s(speaker_scaled_amp_max):VowelLOT 1.004 1.008 30.193 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelNURSE 1.002 1.003 0.342 0.558576
## s(speaker_scaled_amp_max):VowelSTART 1.005 1.009 4.070 0.042796 *
## s(speaker_scaled_amp_max):VowelSTRUT 1.001 1.002 7.153 0.007461 **
## s(speaker_scaled_amp_max):VowelTHOUGHT 3.628 4.539 85.004 < 2e-16 ***
## s(speaker_scaled_amp_max):VowelTRAP 1.271 1.497 7.367 0.007482 **
## s(speaker_scaled_pitch):VowelDRESS 3.349 4.202 5.393 0.000227 ***
## s(speaker_scaled_pitch):VowelFLEECE 1.001 1.003 8.044 0.004520 **
## s(speaker_scaled_pitch):VowelFOOT 1.001 1.002 0.048 0.827860
## s(speaker_scaled_pitch):VowelGOOSE 1.005 1.010 14.754 0.000116 ***
## s(speaker_scaled_pitch):VowelKIT 3.223 4.041 15.382 < 2e-16 ***
## s(speaker_scaled_pitch):VowelLOT 4.270 5.276 9.567 < 2e-16 ***
## s(speaker_scaled_pitch):VowelNURSE 1.006 1.013 1.776 0.180196
## s(speaker_scaled_pitch):VowelSTART 1.001 1.003 5.222 0.022216 *
## s(speaker_scaled_pitch):VowelSTRUT 2.101 2.672 13.202 3.13e-07 ***
## s(speaker_scaled_pitch):VowelTHOUGHT 6.234 7.341 27.210 < 2e-16 ***
## s(speaker_scaled_pitch):VowelTRAP 2.878 3.635 2.562 0.034105 *
## s(Speaker) 134.443 169.000 1505.932 < 2e-16 ***
## s(Vowel,Speaker) 1552.576 1843.000 94.256 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.896 Deviance explained = 89.8%
## fREML = 5.0745e+05 Scale est. = 24994 n = 77901
We inspect the plots for amplitude, articulation rate, and pitch.
all_vowel_plot_f2 <- function(predictions, xvar, xlabel, plot_title) {
vowel_colours_reordered <- c(
FLEECE = "#D89000",
DRESS = "#9590FF",
TRAP = "#FF62BC",
NURSE = "#00BFC4",
GOOSE = "#A3A500",
KIT = "#39B600",
START = "#00B0F6",
STRUT = "#F8766D",
FOOT = "#966432",
LOT = "#00BF7D",
THOUGHT = "#E76BF3"
)
xvar <- enquo(xvar)
predictions %>%
mutate(# Enable plot to distinguish between high and low vowels (for faceting)
height = if_else(Vowel %in% front_vowels, "front", "back or mid")
) %>%
ggplot(
aes(
y = fit,
x = !!xvar,
colour = Vowel
)
) +
geom_ribbon(
aes(
ymin = fit - CI,
ymax = fit + CI,
fill = Vowel
),
alpha = 0.3,
colour = NA
) +
geom_line() +
scale_colour_manual(
values = vowel_colours_reordered
) +
scale_fill_manual(
values = vowel_colours_reordered
) +
labs(
title = plot_title,
x = xlabel,
y = "F2"
)
# facet_wrap(
# facets=vars(height), scales = "free"
# )
}
gamm_preds <- get_predictions(
gamm_fit_f2,
cond = list(
'speaker_scaled_amp_max' = seq(-2.5, 2.5, length.out=200),
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
all_vowel_plot_f2(
gamm_preds,
speaker_scaled_amp_max,
'Scaled amplitude',
'F2 and Amplitude by Vowel'
)
Figure 2.32: F2 and Amplitude by Vowel
thought (the smooth closest to the bottom of the plot) and goose F2 (the brown smooth third from the top, on the left) stick out. This makes sense given our PCA analysis (that is, they are the vowels whose F2 values most strongly pattern against the F1 values in PC1). In general, we see a reduction in F1, although not as large as in the F1 case.
Let’s look at articulation rate:
gamm_preds <- get_predictions(
gamm_fit_f2,
cond = list(
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = seq(-2.5, 2.5, length.out=200),
'speaker_scaled_time' = 0.5,
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
all_vowel_plot_f2(
gamm_preds,
speaker_scaled_art_rate,
'Scaled articulation rate',
'F2 and Articulation Rate by Vowel'
)
Figure 2.33: F2 and Articulation Rate by Vowel
Only dress, fleece, and thought indicated as significant at 0.05. In case of thought, we might have a gradual increase with some surprising ‘wiggles’. dress and fleece F2 seem to decrease as articulation rate increases. This is again consistent with vowel space area reduction.
gamm_preds <- get_predictions(
gamm_fit_f2,
cond = list(
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F',
'speaker_scaled_pitch' = seq(-2.5, 2.5, length.out=200),
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'Vowel' = unique(qb_vowels$Vowel)
)
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; with 200 values ranging from -2.500000 to 2.500000.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
all_vowel_plot_f2(
gamm_preds,
speaker_scaled_pitch,
'Scaled pitch',
'F2 and Pitch by Vowel'
)
Figure 2.34: F2 and Pitch by Vowel
thought F2 is listed in the model summary as significant at 0.05, the plot shows that the smooth deviates from the null hypothesis by wiggling in the middle. It is hard to be confident that this is a real phenomenon.
We now plot the impact of amplitude and articulation rate within the
vowel space as whole. In order to do this, we use the fREML models of
F1 and F2 with the same structure as used in the ML model reported in
the paper.
We first collect the data for amplitude and articulation rate:
# Extends beyond 2.5 as this seems to help the stability of animations below.
plotting_range = seq(-2.6, 2.6, by = 0.1)
amp_gamm_preds_f1 <- get_predictions(
gamm_fit,
cond = list(
# We go slightly beyond the bounds of the model to fix a problem where
# text labels disappear in the animation below.
'speaker_scaled_amp_max' = plotting_range,
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'Vowel' = vowels
),
se = FALSE
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 53 values ranging from -2.600000 to 2.600000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
art_gamm_preds_f1 <- get_predictions(
gamm_fit,
cond = list(
# We go slightly beyond the bounds of the model to fix a problem where
# text labels disappear in the animation below.
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = plotting_range,
'speaker_scaled_time' = 0.5,
'Vowel' = vowels
),
se = FALSE
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 53 values ranging from -2.600000 to 2.600000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
amp_gamm_preds_f2 <- get_predictions(
gamm_fit_f2,
cond = list(
'speaker_scaled_amp_max' = plotting_range,
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = 0,
'speaker_scaled_time' = 0.5,
'Vowel' = vowels
),
se = FALSE
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; set to the value(s): 0.
## * speaker_scaled_amp_max : numeric predictor; with 53 values ranging from -2.600000 to 2.600000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
art_gamm_preds_f2 <- get_predictions(
gamm_fit_f2,
cond = list(
'speaker_scaled_amp_max' = 0,
'participant_gender' = 'F',
'speaker_scaled_pitch' = 0,
'speaker_scaled_art_rate' = plotting_range,
'speaker_scaled_time' = 0.5,
'Vowel' = vowels
),
se = FALSE
)
## Summary:
## * participant_gender : factor; set to the value(s): F.
## * Vowel : factor with 11 values; set to the value(s): DRESS, FLEECE, FOOT, GOOSE, KIT, LOT, NURSE, START, STRUT, THOUGHT, ... (Might be canceled as random effect, check below.)
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.5.
## * speaker_scaled_art_rate : numeric predictor; with 53 values ranging from -2.600000 to 2.600000.
## * speaker_scaled_amp_max : numeric predictor; set to the value(s): 0.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): 0.
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Speaker),s(Vowel,Speaker)
##
gamm_preds_f1 <- bind_rows(
"Articulation Rate" = art_gamm_preds_f1 %>%
select(
Vowel, speaker_scaled_art_rate, fit
) %>%
rename(
scaled_value = speaker_scaled_art_rate
),
"Amplitude" = amp_gamm_preds_f1 %>%
select(
Vowel, speaker_scaled_amp_max, fit
) %>%
rename(
scaled_value = speaker_scaled_amp_max
),
.id = "Variable"
)
gamm_preds_f2 <- bind_rows(
"Articulation Rate" = art_gamm_preds_f2 %>%
select(
Vowel, speaker_scaled_art_rate, fit
) %>%
rename(
scaled_value = speaker_scaled_art_rate
),
"Amplitude" = amp_gamm_preds_f2 %>%
select(
Vowel, speaker_scaled_amp_max, fit
) %>%
rename(
scaled_value = speaker_scaled_amp_max
),
.id = "Variable"
)
gamm_preds <- bind_rows(
"F1" = gamm_preds_f1,
"F2" = gamm_preds_f2,
.id = "Formant"
)
gamm_preds <- gamm_preds %>%
pivot_wider(
names_from = Formant,
values_from = fit
)
first_obs <- gamm_preds %>%
group_by(Vowel, Variable) %>%
slice(which.min(scaled_value))
We now produce a static plot for the paper:
static_plot <- gamm_preds %>%
filter(
between(scaled_value, -2.5, 2.5)
) %>%
ggplot(
aes(
x = F2,
y = F1,
colour = Vowel,
group = Vowel,
label = Vowel,
# frame is introduced here for the interactive plot below. It is
# not necessary for the static plot.
frame = scaled_value
)
) +
geom_path(
arrow = arrow(length = unit(1.5, "mm"), type = "closed"), # Make arrows smaller
show.legend = FALSE
) +
# geom_label(
# fontface = 2,
# size = 2,
# show.legend = FALSE,
# data = first_obs,
# alpha = 0.5
# ) +
geom_point(data=first_obs) +
#label the axes
xlab("F2 (Hz)") +
ylab("F1 (Hz)") +
#reverse the axes to follow conventional vowel plotting
scale_x_reverse(expand = expansion(add=100), position = "top") +
scale_y_reverse(expand = expansion(add=50), position = "right") +
#set the colours. Use of 'rev' means the vowels are labelled roughly from top
#to bottom.
scale_color_manual(values = rev(vowel_colours_with_foot)) +
#add a title
labs(
title = "Vowel Space Effect of Amplitude and Articulation Rate",
subtitle = "By-Speaker Z-Scores Between -2.5 and 2.5"
) +
#set the theme
# theme_bw() +
#make text more visible
theme(
plot.title = element_text(size = 16, hjust = 0, face = "bold"),
plot.subtitle = element_text(size = 16, hjust = 0),
panel.spacing.x = unit(6, "mm"),
strip.text = element_text(size = 10),
# axis.title = element_text(size = 14, face = "bold"),
axis.text.x = element_text(size = 8, face = "bold"),
axis.text.y = element_text(size = 8, face = "bold", angle = 270),
plot.margin = margin(5, 5, 5, 5, "mm"),
# axis.ticks = element_blank(),
# plot.caption = element_text(size = 14, hjust = 0),
# legend.position = "none"
) +
facet_grid(cols = vars(Variable))
ggsave(
here('plots', 'amp-art.png'),
plot = static_plot,
units = "cm",
width = 25,
height = 20
)
static_plot
Figure 2.35: Static plot of effect of amplitude and articulation rate on vowel space as predicted by GAMM models.
We use points rather than labels here because the shifts in articulation rate are so small.
We now produce the same plot as an animation. This will be used in presentations.
animated_plot <- gamm_preds %>%
ggplot(
aes(
x = F2,
y = F1,
colour = Vowel,
group = Vowel,
label = Vowel
)
) +
geom_path(
arrow = arrow(length = unit(1, "mm"), type = "closed"), # Make arrows smaller
show.legend = FALSE
) +
geom_label(
fontface = 2,
size = 2.5,
show.legend = FALSE,
alpha = 0.5
) +
#label the axes
xlab("F2 (Hz)") +
ylab("F1 (Hz)") +
#reverse the axes to follow conventional vowel plotting
scale_x_reverse(expand = expansion(add=200), position = "top") +
scale_y_reverse(expand = expansion(add=50), position = "right") +
#set the colours
scale_color_manual(values = rev(vowel_colours_with_foot)) +
#add a title
labs(
title = "Vowel Space Effect of Amplitude and Articulation Rate",
subtitle = "By-Speaker Z-Scores Between -2.5 and 2.5",
caption = 'Z-scored value: {round(frame_along, 1)}'
) +
#set the theme
# theme_bw() +
#make text more visible
theme(
plot.title = element_text(size = 16, hjust = 0, face = "bold"),
plot.subtitle = element_text(size = 16, hjust = 0),
axis.title = element_text(size = 14, face = "bold"),
axis.text.x = element_text(size = 14, face = "bold"),
axis.text.y = element_text(size = 14, face = "bold", angle = 270),
axis.ticks = element_blank(),
plot.caption = element_text(size = 14, hjust = 0),
# legend.position = "none"
) +
facet_grid(cols = vars(Variable)) +
transition_reveal(along = scaled_value, range = c(-2.5, 2.5))
animated_plot <- animate(
animated_plot,
start_pause = 10,
end_pause = 10
)
animated_plot
Figure 2.36: Animated plot of effect of amplitude and articulation rate on vowel space as predicted by GAMM models.
Finally, we produce an interactive plot using plotly.
interactive_plot <- gamm_preds %>%
filter(
between(scaled_value, -2.5, 2.5)
) %>%
ggplot(
aes(
x = F2,
y = F1,
colour = Vowel,
group = Vowel,
label = Vowel,
# frame is introduced here for the interactive plot below. It is
# not necessary for the static plot.
frame = scaled_value
)
) +
geom_text(
fontface = 2,
size = 3,
show.legend = FALSE,
alpha = 1
) +
#label the axes
xlab("F2 (Hz)") +
ylab("F1 (Hz)") +
#reverse the axes to follow conventional vowel plotting
scale_x_reverse(expand = expansion(add=200), position = "top") +
scale_y_reverse(expand = expansion(add=50), position = "right") +
#set the colours. Use of 'rev' means the vowels are labelled roughly from top
#to bottom.
scale_color_manual(values = rev(vowel_colours_with_foot)) +
#add a title
labs(
title = "Vowel Space Effect of Amplitude and Articulation Rate",
subtitle = "By-Speaker Z-Scores Between -2.5 and 2.5"
) +
#set the theme
# theme_bw() +
#make text more visible
theme(
plot.title = element_text(size = 16, hjust = 0, face = "bold"),
plot.subtitle = element_text(size = 16, hjust = 0),
axis.title = element_text(size = 14, face = "bold"),
axis.text.x = element_text(size = 14, face = "bold"),
axis.text.y = element_text(size = 14, face = "bold", angle = 270),
axis.ticks = element_blank(),
plot.caption = element_text(size = 14, hjust = 0),
legend.position = "none"
) +
facet_grid(cols = vars(Variable))
ggplotly(interactive_plot)
Figure 2.37: Interactive plot of effect of amplitude and articulation rate on vowel space as predicted by GAMM models.
The takeaway from these plots is simple: the widely controlled for effect of articulation rate on the vowel space is much smaller than the effect of maximum amplitude.
We now turn to whether variation in relative amplitude, perhaps through changes to articulatory setting, is being used by speakers to indicate topical structure in the QuakeBox monologues.
We are interested in whether there is a statistically significant effect of coming to the end of a topic and a change in amplitude. The most sophisticated way to investigate this is to consider the trajectory of amplitude over time with a flexible nonlinear model like a GAMM. However, this level of sophistication introduces all sorts of complexities. It is worth starting with a simpler model and seeing what we find. One very simple model assigns an amplitude value to the beginning, middle, and end of each topic chunk. We will then see whether there is any systematic difference between these values in our corpus.
We begin by dividing the topics into three even (by time through the topic) chunks.
qb_vowels <- qb_vowels %>%
group_by(Speaker, topic_no) %>%
mutate(
topic_part = cut(
topic_time_scaled,
breaks = c(-0.1, 0.33, 0.66, 1.1),
labels = c("start", "middle", "end")
)
) %>%
group_by(Speaker, topic_no, topic_part) %>%
mutate(
topic_part_n = n()
) %>%
ungroup()
We check that there is an even amount of tokens in the beginning, middle, and end chunks (Figure 3.4). Visually, the boxes look reasonably similar.
qb_vowels %>%
ggplot(
aes(
x = topic_part,
y = topic_part_n
)
) +
geom_boxplot() +
labs(
title = "Topic Parts by Token Count",
x = "Topic part",
y = "Token count"
)
Figure 3.4: Topic parts by token count.
Given the code above, it is not possible to have a topic_part_n score of 0,
but it may be that there are some empty parts. We create a new dataframe which
will have a zero value for any topic part without tokens.
qb_parts <- qb_vowels %>%
group_by(Speaker, topic_no, topic_part, .drop=FALSE) %>%
summarise(
topic_part_n = first(topic_part_n),
topic_part_n = if_else(is.na(topic_part_n), 0L, topic_part_n) # Replace NAs with 0
)
We now look at a historgram of the number of tokens within topic parts (Figure 3.5).
qb_parts %>%
ggplot(
aes(
x = topic_part_n
)
) +
geom_histogram(bins = 100) +
labs(
title = "Distribution of Token Counts Within Each Topic Part.",
y = "Count of topic parts.",
x = "Count of tokens within topic part."
)
Figure 3.5: Distribution of token counts within each topic part.
There are not many parts without a token.
Our linear model will fit a mean to each of our topic parts for each speaker and topic. A simple way to ensure enough data for this is to insist on 5 tokens, as we did for our initial speaker filtering.
topics_to_filter <- qb_vowels %>%
filter(
topic_part_n < 5
) %>%
select(Speaker, topic_no) %>%
unique()
speaker_topics_to_filter <- topics_to_filter %>%
mutate(
speaker_topic = str_c(Speaker, "_", topic_no)
) %>%
pull(speaker_topic)
qb_filtered <- qb_vowels %>%
mutate(
speaker_topic = str_c(Speaker, "_", topic_no)
) %>%
filter(
!speaker_topic %in% speaker_topics_to_filter
)
We lose 63 topics from the data set with this filtering (from 617).
For convenience, we keep only those variables that we want for modelling.
qb_glmm_data <- qb_filtered %>%
select(
Speaker, participant_gender, participant_age_category,
Vowel, speaker_scaled_art_rate, speaker_topic, topic_part,
topic_time_scaled, speaker_scaled_time, speaker_scaled_amp_max,
speaker_scaled_pitch, topic_no, topic, time, speaker_length
)
We now look to see if there is any evidence that the parts of topics are systematically different in terms of amplitude.
Before formally modelling, we produce a boxplot of amplitude at the start middle and end of topics. (Figure 3.6).
qb_glmm_data %>%
ggplot(
aes(
x = topic_part,
y = speaker_scaled_amp_max
)
) +
geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), alpha=0.2) +
geom_jitter(alpha=0.02) +
labs(
title = "Start, Middle, and End Mean Amplitudes",
y = "Scaled maximum amplitude",
x = "Part of topic"
)
Figure 3.6: Topic start, middle, and end amplitude
We do seem to have a small amplitude reduction at the end of a topic. Since it is hard to read the precise values off Figure 3.6, we have a look at the means:
qb_glmm_data %>%
group_by(topic_part) %>%
summarise(speaker_scaled_amp_max = mean(speaker_scaled_amp_max, na.rm=TRUE))
The start and middle sit slightly above the average amplitude, while the end sits below. However, we know that amplitude tends to be below the mean at the end of monologues and this effect may hide the relative shifts in amplitude within topical chunks of the monologue.
We fit a model of amplitude with topic part, scaled time through the speaker’s
monologue, and pitch scaled at the speaker level as predictors. We allow each
speaker to have their own speaker_scaled_time slope to capture the overall
trend in amplitude in their monologue. We do not give each speaker a random
intercept, since the data is scaled for each speaker. However, we centre the
time values around 0, so that the intercept is in the centre of the monologue
and changes in slope can more effectively capture changes over the course of the
monologue. We also treat each topical chunk independently, allowing the
beginning, middle, and end to be fit separately for each topical section (as a
random effect). We could treat topic part as an ordered factor, but I have stuck
with a simpler model which allows each part to vary independently.
NB: Change eval=FALSE to eval=TRUE to fit this model yourself.
glmm_fit <- lmer(
speaker_scaled_amp_max ~
topic_part +
speaker_scaled_time +
speaker_scaled_pitch +
(0+speaker_scaled_time|Speaker) +
(0+topic_part|speaker_topic),
data=qb_glmm_data %>%
mutate(
speaker_scaled_time = speaker_scaled_time - 0.5
))
# This model needs some help converging. We use the strategy of restarting
# from current parameters
pars <- getME(glmm_fit, "theta")
glmm_fit <- update(glmm_fit, start=pars)
write_rds(glmm_fit, here('models', 'glmm_fit_pitch.rds'))
We look at the summary.
glmm_fit <- read_rds(here('models', 'glmm_fit_pitch.rds'))
summary(glmm_fit)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## speaker_scaled_amp_max ~ topic_part + speaker_scaled_time + speaker_scaled_pitch +
## (0 + speaker_scaled_time | Speaker) + (0 + topic_part | speaker_topic)
## Data: qb_glmm_data %>% mutate(speaker_scaled_time = speaker_scaled_time -
## 0.5)
##
## REML criterion at convergence: 171076
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5987 -0.6476 -0.0052 0.6522 4.3801
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## speaker_topic topic_partstart 0.04660 0.2159
## topic_partmiddle 0.03713 0.1927 0.44
## topic_partend 0.04835 0.2199 0.21 0.51
## Speaker speaker_scaled_time 0.21681 0.4656
## Residual 0.55657 0.7460
## Number of obs: 74978, groups: speaker_topic, 554; Speaker, 153
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.059289 0.011389 5.206
## topic_partmiddle -0.010110 0.012815 -0.789
## topic_partend -0.071852 0.015365 -4.676
## speaker_scaled_time -0.203193 0.044822 -4.533
## speaker_scaled_pitch 0.541694 0.003559 152.190
##
## Correlation of Fixed Effects:
## (Intr) tpc_prtm tpc_prtn spkr_scld_t
## tpc_prtmddl -0.636
## topic_prtnd -0.655 0.626
## spkr_scld_t 0.104 -0.103 -0.175
## spkr_scld_p 0.005 0.006 0.017 0.012
This model provides some evidence for a reduction in amplitude at the end of each topical chunk which is not explained by the reduction in amplitude over the course of the monologue as a whole or by variation in pitch. We don’t get an explicit p-value, but our t values for a higher than average amplitude at the start of the monologue and a lower than average amplitude at the end of a monologue have t-values much greater than two, which we take to be surprising.3
The magnitude of the drop at the end of a topical span is somewhat larger than
the increase in amplitude at the start of a topical span. The reduction in
amplitude over the course of an entire monologue can be easily read off the
speaker_scaled_time coefficient, which says an increase in one in scaled time
is associated with a reduction in amplitude of -0.24 (scaled amplitude).
Obviously, an increase in monologue length of one represents going through the
entire monologue.
We have a look at the QQ-plot:
qqnorm(resid(glmm_fit))
qqline(resid(glmm_fit), col=2)
Figure 3.7: QQ Plot for GLMM model
Figure 3.7 shows that our residuals are behaving (sufficiently) normally.
Another way to try to handle the known reduction in amplitude over the course of a monologue would be to remove all end-of-monologue topics. We filter out all monologues which end at the same time that the monologue ends.
start_time_filter <- qb_glmm_data %>%
group_by(Speaker, topic_no) %>%
summarise(
end_time = max(time),
speaker_length = first(speaker_length)
) %>%
mutate(
speaker_topic = str_c(Speaker, "_", topic_no)
) %>%
filter(
end_time == speaker_length
) %>%
pull(speaker_topic)
glmm_fit_trimmed <- lmer(
speaker_scaled_amp_max ~
topic_part +
speaker_scaled_time +
speaker_scaled_pitch +
(0+speaker_scaled_time|Speaker) +
(0+topic_part|speaker_topic),
data=qb_glmm_data %>%
filter(!speaker_topic %in% start_time_filter) %>%
mutate(
speaker_scaled_time = speaker_scaled_time - 0.5
))
write_rds(glmm_fit_trimmed, here('models', 'glmm_fit_trimmed.rds'))
glmm_fit_timmed <- read_rds(here('models', 'glmm_fit_trimmed.rds'))
summary(glmm_fit_timmed)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## speaker_scaled_amp_max ~ topic_part + speaker_scaled_time + speaker_scaled_pitch +
## (0 + speaker_scaled_time | Speaker) + (0 + topic_part | speaker_topic)
## Data: qb_glmm_data %>% filter(!speaker_topic %in% start_time_filter) %>%
## mutate(speaker_scaled_time = speaker_scaled_time - 0.5)
##
## REML criterion at convergence: 144408.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4849 -0.6459 -0.0060 0.6519 4.4243
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## speaker_topic topic_partstart 0.04668 0.2161
## topic_partmiddle 0.03771 0.1942 0.41
## topic_partend 0.04592 0.2143 0.18 0.47
## Speaker speaker_scaled_time 0.21471 0.4634
## Residual 0.55305 0.7437
## Number of obs: 63457, groups: speaker_topic, 474; Speaker, 143
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.076129 0.012679 6.004
## topic_partmiddle -0.019891 0.014087 -1.412
## topic_partend -0.080855 0.016637 -4.860
## speaker_scaled_time -0.186784 0.047877 -3.901
## speaker_scaled_pitch 0.555734 0.003874 143.434
##
## Correlation of Fixed Effects:
## (Intr) tpc_prtm tpc_prtn spkr_scld_t
## tpc_prtmddl -0.640
## topic_prtnd -0.684 0.618
## spkr_scld_t 0.177 -0.108 -0.185
## spkr_scld_p 0.007 0.004 0.016 0.012
The results of the two models seem roughly equivalent. Since the difference between the two isn’t great, we’ll just visualise the first model.
We set up a data frame for the data which we want to draw predictions from. We
set speaker_time_scaled to 0, which now represents the middle of the monologue.
That is, we are asking the model to assume that the speaker is half way through
their monologue.
new_data <- tibble(
speaker_topic = rep(unique(qb_glmm_data$speaker_topic), times=3),
topic_part = rep(
c('start', 'middle', 'end'),
each=length(unique(qb_glmm_data$speaker_topic))
),
# Assuming we are in the middle of the monologue
speaker_scaled_time = 0,
# and at average pitch
speaker_scaled_pitch = 0
)
new_data <- new_data %>%
filter(
!is.na(speaker_topic)
) %>%
mutate(
Speaker = str_extract(speaker_topic, 'QB_NZ_[MF]_[0-9]+')
)
new_data <- new_data %>%
mutate(
prediction = predict(glmm_fit, newdata=new_data)
)
We plot the model estimates for each topic in the dataset in Figure 3.8.
new_data %>%
mutate(
topic_part = factor(topic_part, levels = c('start', 'middle', 'end'))
) %>%
ggplot(
aes(
x = topic_part,
y = prediction
)
) +
geom_violin(
draw_quantiles = c(0.25, 0.5, 0.75),
alpha = 0.5
) +
geom_jitter(
alpha = 0.2
) +
labs(
title = "Model Predictions for Each Topic",
x = "Part",
y = "Predicted amplutide"
)
Figure 3.8: Model predictions for each topic at mid point of monologue.
We see a clear difference in the means of the different parts. We will generate confidence intervals for our mean estimates after carrying out a check using simulated topics below.
An interesting question which this data can address: do people talk relatively quietly when discussing different topics? As noted above, this is not our main focus, but we can look at relevant data in passing.
qb_glmm_data %>%
group_by(Speaker) %>%
ggplot(
aes(
x = topic,
y = speaker_scaled_amp_max
)
) +
geom_boxplot() +
theme(axis.text.x = element_text(angle=90)) +
labs(
title = "Scaled Amplitude by Topic.",
x = "Topic",
y = "Scaled amplitude"
)
Figure 3.9: Scaled amplitude by topic.
Figure 3.9 suggests lower amplitude when discussing the December earthquake (only 17 speakers) and the June earthquake (50 speakers), and a higher amplitude when discussing thoughts for the future.
We fit a model to address the same question.
glmm_fit_topic <- lmer(
speaker_scaled_amp_max ~
topic +
speaker_scaled_time +
speaker_scaled_pitch +
(0+speaker_scaled_time|Speaker) +
(1|speaker_topic),
data=qb_glmm_data %>%
mutate(
speaker_scaled_time = speaker_scaled_time - 0.5
)
)
summary(glmm_fit_topic)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## speaker_scaled_amp_max ~ topic + speaker_scaled_time + speaker_scaled_pitch +
## (0 + speaker_scaled_time | Speaker) + (1 | speaker_topic)
## Data: qb_glmm_data %>% mutate(speaker_scaled_time = speaker_scaled_time -
## 0.5)
##
## REML criterion at convergence: 172012.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8041 -0.6508 -0.0032 0.6570 4.5562
##
## Random effects:
## Groups Name Variance Std.Dev.
## speaker_topic (Intercept) 0.02904 0.1704
## Speaker speaker_scaled_time 0.27469 0.5241
## Residual 0.57034 0.7552
## Number of obs: 74978, groups: speaker_topic, 554; Speaker, 153
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.013943 0.017525 0.796
## topic{december 23rd} 0.015257 0.057408 0.266
## topic{february} 0.064615 0.023259 2.778
## topic{government response} 0.001543 0.058470 0.026
## topic{housing and insurance} -0.016842 0.041389 -0.407
## topic{june} 0.096088 0.032719 2.937
## topic{personal background} -0.062548 0.044957 -1.391
## topic{september} -0.028893 0.029059 -0.994
## topic{thoughts for the future} 0.137887 0.090331 1.526
## speaker_scaled_time -0.264335 0.047201 -5.600
## speaker_scaled_pitch 0.542476 0.003564 152.228
##
## Correlation of Fixed Effects:
## (Intr) tp{23} tpc{f} tpc{r} tp{ai} tpc{j} tpc{b} tpc{s} t{ftf}
## tpc{dc23rd} -0.256
## topc{fbrry} -0.768 0.187
## tpc{grspns} -0.282 0.078 0.211
## tpc{ainsrn} -0.391 0.112 0.295 0.118
## topic{june} -0.481 0.159 0.363 0.143 0.201
## tpc{bckgrn} -0.426 0.103 0.330 0.113 0.152 0.184
## tpc{sptmbr} -0.703 0.140 0.547 0.183 0.251 0.281 0.324
## tpc{tftftr} -0.174 0.052 0.125 0.046 0.069 0.090 0.059 0.098
## spkr_scld_t -0.132 -0.017 0.117 0.018 0.020 -0.011 0.095 0.207 -0.010
## spkr_scld_p -0.001 0.005 0.008 0.006 0.006 -0.003 0.012 0.013 0.003
## spkr_scld_t
## tpc{dc23rd}
## topc{fbrry}
## tpc{grspns}
## tpc{ainsrn}
## topic{june}
## tpc{bckgrn}
## tpc{sptmbr}
## tpc{tftftr}
## spkr_scld_t
## spkr_scld_p 0.018
This gives no clear evidence for or against any significant effect here. Exploring whether these different specific topics have any discernible phonetic features is left for future work.
Another way to proceed is to fit GAMM models over topic windows with a
speaker_scaled_time control variable. We’ll continue with the data filtered
for the three part GLMM approach as data gaps can lead to poor behaviour from
GAMMs.
Let’s look at an uncontrolled version produced using ggplot’s geom_smooth.
qb_filtered %>%
group_by(Speaker) %>%
ggplot(
aes(
x = topic_time_scaled,
y = speaker_scaled_amp_max
)
) +
geom_smooth() +
labs(
title = "Scaled amplitude over course of topic.",
y = "Scaled amplitude",
x = "Scaled topic time"
)
Figure 3.10: Smooth of scaled amp by topic time.
Figure 3.10 shows a very similar pattern to Figure 2.13. That is, we get a similar drop over the course of a topic as we get over the course of the monologue as a whole.
We can, again as a side track, ask if the pattern is different for different topics.
qb_filtered %>%
group_by(Speaker) %>%
ggplot(
aes(
x = topic_time_scaled,
y = speaker_scaled_amp_max
)
) +
geom_smooth() +
facet_wrap(vars(topic)) +
labs(
title = "Scaled amplitude over course of topic (specific).",
y = "Scaled amplitude",
x = "Scaled topic time"
)
Figure 3.11: Smooth of scaled amp by topic time, by specific topic.
The relative size of the confidence intervals here mostly tracks different
amounts of data available for each topic. Those with lots of data seem to have
gradual slow declines in amplitude. Certainly, there are none with clear and
systematic increases over the course of the topic. The drop off seems to start
late, rather than just being a straight line. This is clearest in the case of
the Feb and Sept earthquakes. Perhaps {aftermath} is a counterexample to this.
Again, this is left open for future research.
It is worth noting that here we are particularly at risk of of misinterpreting the reduction in amplitude over the entire monologue as an effect of topic. Some topics, like the February and September earthquakes, take up large portions of speaker monologues so we are likely to see any overall reduction come through within these topical sections even if the reduction has nothing to do with the topics themselves.
We now model using a GAMM approach. We allow random smooths for each topic in the corpus.
We turn speaker_topic into an R factor for the purpose of fitting the GAMM.
qb_gamm_data <- qb_glmm_data %>%
mutate(
speaker_topic = as.factor(speaker_topic)
)
gamm_topic <- bam(
speaker_scaled_amp_max ~
s(topic_time_scaled) +
s(speaker_scaled_time) +
s(speaker_scaled_pitch) +
s(topic_time_scaled, speaker_topic, bs = "fs", k=5, m=1) +
s(speaker_scaled_time, Speaker, bs = 'fs', k=5, m=1),
data = qb_gamm_data,
method = 'fREML',
discrete = TRUE,
nthreads = 8
)
write_rds(gamm_topic, here('models', 'gamm_topic.rds'))
gamm_topic_summary <- summary(gamm_topic)
write_rds(gamm_topic_summary, here('models', 'gamm_topic_summary.rds'))
gamm_topic <- read_rds(here('models', 'gamm_topic.rds'))
gamm_topic_summary <- read_rds(here('models', 'gamm_topic_summary.rds'))
gamm_topic_summary
##
## Family: gaussian
## Link function: identity
##
## Formula:
## speaker_scaled_amp_max ~ s(topic_time_scaled) + s(speaker_scaled_time) +
## s(speaker_scaled_pitch) + s(topic_time_scaled, speaker_topic,
## bs = "fs", k = 5, m = 1) + s(speaker_scaled_time, Speaker,
## bs = "fs", k = 5, m = 1)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.08458 0.00715 11.83 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(topic_time_scaled) 6.687 7.783 9.878 <2e-16 ***
## s(speaker_scaled_time) 6.733 7.516 11.386 <2e-16 ***
## s(speaker_scaled_pitch) 8.552 8.932 2923.574 <2e-16 ***
## s(topic_time_scaled,speaker_topic) 1155.415 2769.000 1.102 <2e-16 ***
## s(speaker_scaled_time,Speaker) 259.722 764.000 1.111 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.336 Deviance explained = 34.8%
## fREML = 84021 Scale est. = 0.53007 n = 74978
We check the model diagnostics.
gam.check(gamm_topic)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] 1.220064e-10 -9.290892e-10 1.643401e-06 -8.127131e-09 -5.865115e-08
## [6] 9.081560e-09 -3.316299e-04 -1.820110e-06
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## 1.712436e+00 8.566236e-02 3.147857e-04 2.020869e-01 2.564236e-02
## 8.566236e-02 1.966183e+00 -1.869969e-03 1.378672e-01 -2.612888e-01
## 3.147857e-04 -1.869969e-03 3.680827e+00 -1.706017e-02 -1.324081e-01
## 2.020869e-01 1.378672e-01 -1.706017e-02 1.993749e+02 1.262124e+00
## 2.564236e-02 -2.612888e-01 -1.324081e-01 1.262124e+00 5.671660e+01
## 5.667508e-02 7.754806e-02 1.921780e-02 9.965837e+00 1.916302e+01
## -1.067414e-07 3.700524e-07 -4.468810e-06 1.423235e-05 2.796439e-04
## d -2.843337e+00 -2.866275e+00 -3.776014e+00 -4.535389e+02 -1.241687e+02
## [,6] [,7] [,8]
## 5.667508e-02 -1.067414e-07 -2.843337e+00
## 7.754806e-02 3.700524e-07 -2.866275e+00
## 1.921780e-02 -4.468810e-06 -3.776014e+00
## 9.965837e+00 1.423235e-05 -4.535389e+02
## 1.916302e+01 2.796439e-04 -1.241687e+02
## 7.601142e+01 2.502130e-05 -1.298600e+02
## 2.502130e-05 3.316399e-04 -8.256806e-04
## d -1.298600e+02 -8.256806e-04 3.748700e+04
##
## Model rank = 3563 / 3563
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(topic_time_scaled) 9.00 6.69 1.01 0.62
## s(speaker_scaled_time) 9.00 6.73 0.98 0.12
## s(speaker_scaled_pitch) 9.00 8.55 1.01 0.70
## s(topic_time_scaled,speaker_topic) 2770.00 1155.42 1.01 0.64
## s(speaker_scaled_time,Speaker) 765.00 259.72 0.98 0.13
It looks like our default k values are OK. All other diagnostics are fine.
Let’s have a look at the smooths:
plot_smooth(gamm_topic, view="topic_time_scaled")
## Summary:
## * topic_time_scaled : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.486097114063216.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): -0.067742628346803.
## * speaker_topic : factor; set to the value(s): QB_NZ_F_408_2. (Might be canceled as random effect, check below.)
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(topic_time_scaled,speaker_topic),s(speaker_scaled_time,Speaker)
##
Figure 3.12: Amplitude by time through topic
We see a strong uptick at the start and a strong downward movement at the end of topics. Here the model is assuming that we are in the middle of a monologue.
We now look at amplitude through the monologue as predicted by this model.
plot_smooth(gamm_topic, view="speaker_scaled_time")
## Summary:
## * topic_time_scaled : numeric predictor; set to the value(s): 0.497835510463644.
## * speaker_scaled_time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): -0.067742628346803.
## * speaker_topic : factor; set to the value(s): QB_NZ_F_408_2. (Might be canceled as random effect, check below.)
## * Speaker : factor; set to the value(s): QB_NZ_M_291. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(topic_time_scaled,speaker_topic),s(speaker_scaled_time,Speaker)
##
Figure 3.13: Amplitude by time through Monologue
Here the model assumes we are midway through a topic, and tells us that we drop below the mean value at around 70% of the way through a monologue and then drop to around -0.2 by the end.
Since behaviour of smooths can be affected by data gaps, we look to see that there is a similar amount of data in each part of our topical chunks.
qb_gamm_data %>%
mutate(
topic_parts = cut(topic_time_scaled, breaks=10)
) %>%
group_by(topic_parts) %>%
summarise(
n = n()
)
These values are reasonably consistent with one another.
These results are consistent with those we obtained from our linear mixed models above.
We can see how the different topic smooths look for a particular speaker as follows:
plot_smooth(
gamm_topic,
view="topic_time_scaled",
plot_all="speaker_topic",
cond=list(
speaker_topic = c(
"QB_NZ_F_131_1", "QB_NZ_F_131_2",
"QB_NZ_F_131_3", "QB_NZ_F_131_6",
"QB_NZ_F_131_7"
),
Speaker = c("QB_NZ_F_131")
),
rm.ranef=FALSE)
## Summary:
## * topic_time_scaled : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.486097114063216.
## * speaker_scaled_pitch : numeric predictor; set to the value(s): -0.067742628346803.
## * speaker_topic : factor; set to the value(s): QB_NZ_F_131_1, QB_NZ_F_131_2, QB_NZ_F_131_3, QB_NZ_F_131_6, QB_NZ_F_131_7.
## * Speaker : factor; set to the value(s): QB_NZ_F_131.
All of these topics start near or above the mean and end below the mean, with varying behaviour in the middle.
One way to test whether the drop in amplitude at the end of topical chunks is a real feature of the data or merely explained by the fact that amplitude has a tendency to drop over the monologue, is to create fake topics. That is, we will cut up the monologue into ‘chunks’ with no regard to what the speaker is talking about.
Let’s look at the lengths for each topic in our dataset. We want our fake topical selections to match the real ones in terms of length.
topic_lengths <- qb_vowels %>%
group_by(Speaker, topic_no) %>%
summarise(
topic_length = first(topic_length)
) %>%
filter(
!is.na(topic_no)
)
topic_lengths
Our aim here is to avoid a spurious effect of topic part being generated by
the decrease in amplitude over the course of an entire monologue. The
extent of this effect might be affected if our fake topics are, on the whole
shorter or longer than the real topics. One way to ensure we have the same
distribution of speaker topic lengths is to take random chunks from a
speaker monologue which match the length of their actual topics. So, in the
case of QB_NZ_F_131 we have 7 topic lengths available. Our first fake chunk
would take a random continuous section of the monologue of 731.2 seconds, our
second would be a random continuous section of the monologue of length 96.4.
To do this, we take the lengths of each speaker’s topics and then replace the
real topics from the speaker with a random span of the same length from the
speaker’s monologue. In the following code block the chunks dataframe takes
information about the speakers full monologue length and the length of
each topic, it them sets a random start point for each of the new
chunks (ensuring that there is enough time left in the monologue). The function
collect_chunk is then used to collect the data for the defined spans.
chunks <- qb_vowels %>%
group_by(Speaker, topic_no) %>%
summarise(
topic_length = first(topic_length),
speaker_length = first(speaker_length)
) %>%
filter(
!is.na(topic_no)
) %>%
mutate(
chunk_start = map2_int(
topic_length,
speaker_length,
~ sample(seq(0, round(.y - .x)), 1)
),
chunk_end = chunk_start + topic_length
)
collect_chunk <- function(speaker, chunk_start, chunk_end, vowel_data) {
out_df <- vowel_data %>%
filter(
Speaker == speaker,
time > chunk_start,
time < chunk_end
) %>%
select(
c(
time, art_rate, Vowel, pitch, speaker_scaled_time,
speaker_scaled_amp_max, topic_no, intensity_max,
speaker_scaled_art_rate, speaker_scaled_pitch
)
)
}
chunks <- chunks %>%
mutate(
chunk_df = pmap(
list(Speaker, chunk_start, chunk_end),
~ collect_chunk(..1, ..2, ..3, qb_vowels)
)
) %>%
rename(
chunk = topic_no # We may need to compare overlap of real and fake topics.
) %>%
unnest(chunk_df)
The following blocks repeat the filtering steps carried out above.
We define the start, middle and end and determine how many tokens are in each.
chunks <- chunks %>%
group_by(Speaker, chunk) %>%
mutate(
chunk_time_scaled = rescale(time, to=c(0,1)),
chunk_part = cut(
chunk_time_scaled,
breaks = c(-0.1, 0.33, 0.66, 1.1),
labels = c("start", "middle", "end")
)
) %>%
group_by(Speaker, chunk, chunk_part) %>%
mutate(
chunk_part_n = n()
) %>%
ungroup()
We have a look at how many tokens we have in each chunk part.
chunk_counts <- chunks %>%
group_by(Speaker, chunk, chunk_part, .drop=FALSE) %>%
summarise(
chunk_part_n = first(chunk_part_n),
chunk_part_n = if_else(
is.na(chunk_part_n), 0L, chunk_part_n) # Replace NAs with 0
)
chunk_counts %>%
ggplot(
aes(
x = chunk_part_n
)
) +
geom_histogram(bins = 100) +
labs(
title = "Distribution of Token Counts Within Each Chunk Part.",
y = "Count of chunk parts.",
x = "Count of tokens within chunks part."
)
Figure 3.14: Number of tokens in each fake topic part.
Looks very similar to 3.5. This should not be surprising!
We now remove any chunk with less than five tokens.
chunks_to_filter <- chunks %>%
filter(
chunk_part_n < 5
) %>%
select(Speaker, chunk) %>%
unique()
speaker_chunks_to_filter <- chunks_to_filter %>%
mutate(
speaker_chunk = str_c(Speaker, "_", chunk)
) %>%
pull(speaker_chunk)
chunks_filtered <- chunks %>%
mutate(
speaker_chunk = str_c(Speaker, "_", chunk)
) %>%
filter(
!speaker_chunk %in% speaker_chunks_to_filter
)
We now look at the difference between the lengths of the real and fake topics again (Figure 3.15).
real_topic_lengths <- qb_filtered %>%
group_by(Speaker, topic_no) %>%
summarise(
topic_length = first(topic_length)
)
fake_topic_lengths <- chunks_filtered %>%
group_by(Speaker, chunk) %>%
summarise(
topic_length = first(topic_length) # chunk_length
) %>%
mutate(
Speaker = as.factor(Speaker)
)
bind_rows(
'Real' = real_topic_lengths,
'Fake' = fake_topic_lengths,
.id = "Source"
) %>%
ggplot(
aes(
x = topic_length,
colour = Source
)
) +
geom_freqpoly(bins=100)
Figure 3.15: Length distribution of real and fake topics.
The only purpose of this plot is to ensure that we have basically the same length distribution. The small differences are because different topics/chunks will be filtered out by our filtering rules (ensuring each part has at least 5 tokens).
We now copy the GLMM model specification used on the real topic data.
We also save the current ‘chunks’ and reload them at the next stage. If new chunks are generated above, then different speakers may be filtered out and it will be impossible to use the visualisation code later in this document. Saving the ‘chunks’ at the same time we save the model allows us to ensure that the visualisations below will work. This is unfortunate for reproducibility, but setting a seed at the start of the document does not ensure that the same ‘chunks’ will be generated in interactive use of this code.
glmm_fit_2 <- lmer(
speaker_scaled_amp_max ~
chunk_part +
speaker_scaled_time +
speaker_scaled_pitch +
(0+speaker_scaled_time|Speaker) +
(0+chunk_part|speaker_chunk),
data=chunks_filtered %>%
mutate(
speaker_scaled_time = speaker_scaled_time - 0.5
))
# This model needs some help converging. We use the strategy of restarting
# from current parameters
pars <- getME(glmm_fit_2, "theta")
glmm_fit_2 <- update(glmm_fit_2, start=pars)
write_rds(glmm_fit_2, here('models', 'glmm_pitch_fit_chunks.rds'))
write_rds(chunks_filtered, here('processed_data', 'chunks_filtered.rds'))
glmm_fit_2 <- read_rds(here('models', 'glmm_pitch_fit_chunks.rds'))
chunks_filtered <- read_rds(here('processed_data', 'chunks_filtered.rds'))
summary(glmm_fit_2)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## speaker_scaled_amp_max ~ chunk_part + speaker_scaled_time + speaker_scaled_pitch +
## (0 + speaker_scaled_time | Speaker) + (0 + chunk_part | speaker_chunk)
## Data:
## chunks_filtered %>% mutate(speaker_scaled_time = speaker_scaled_time -
## 0.5)
##
## REML criterion at convergence: 164196.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3559 -0.6470 -0.0008 0.6523 4.3698
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## speaker_chunk chunk_partstart 0.04387 0.2094
## chunk_partmiddle 0.03568 0.1889 0.41
## chunk_partend 0.02586 0.1608 0.02 0.21
## Speaker speaker_scaled_time 0.32373 0.5690
## Residual 0.55997 0.7483
## Number of obs: 71846, groups: speaker_chunk, 519; Speaker, 153
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.069081 0.011944 5.784
## chunk_partmiddle 0.002111 0.013289 0.159
## chunk_partend -0.011741 0.015773 -0.744
## speaker_scaled_time -0.158448 0.054393 -2.913
## speaker_scaled_pitch 0.549359 0.003656 150.260
##
## Correlation of Fixed Effects:
## (Intr) chnk_prtm chnk_prtn spkr_scld_t
## chnk_prtmdd -0.640
## chunk_prtnd -0.771 0.586
## spkr_scld_t 0.153 -0.122 -0.207
## spkr_scld_p 0.019 -0.005 0.001 0.005
The ‘chunks’ certainly provide weaker evidence of a reduction in amplitude with fake topics. This may simply be due to the reduction in amount of data going in to the model (our chunking process is, in some ways, a filtering process). The t-value for the intercept is surprisingly high! Perhaps the model is capturing the pattern of reduction in amplitude by pushing the ‘start’ of the chunks up in amplitude rather than pushing the ‘end’ down. This is not great for our claim that speakers drop amplitude at the end of topical sections. It is not direct counter evidence though. We will find a way to a more principled test below.
We now produce a plot equivalent to Figure 3.8.
new_data_fake <- tibble(
speaker_chunk = rep(unique(chunks_filtered$speaker_chunk), times=3),
chunk_part = rep(
c('start', 'middle', 'end'),
each=length(unique(chunks_filtered$speaker_chunk))
),
# Assuming we are in the middle of the monologue
speaker_scaled_time = 0,
speaker_scaled_pitch = 0
)
new_data_fake <- new_data_fake %>%
filter(
!is.na(speaker_chunk)
) %>%
mutate(
Speaker = str_extract(speaker_chunk, 'QB_NZ_[MF]_[0-9]+')
)
new_data_fake <- new_data_fake %>%
mutate(
prediction = predict(glmm_fit_2, newdata=new_data_fake)
)
new_data_fake %>%
mutate(
# Make sure chunks are in correct order in plot.
chunk_part = factor(chunk_part, levels = c('start', 'middle', 'end'))
) %>%
ggplot(
aes(
x = chunk_part,
y = prediction
)
) +
geom_violin(
draw_quantiles = c(0.25, 0.5, 0.75),
alpha = 0.5
) +
geom_jitter(
alpha = 0.1
) +
labs(
title = "Model Predictions for Each Chunk (Fake Topic) with Position in Monologue Controlled.",
x = "Part",
y = "Predicted amplutide"
)
Figure 3.16: Model predictions for each topic at mid point of monologue.
All are assumed to be above 0 in scaled amplitude in this model, with a slightly increased value for the start.
To make this even clearer, we plot the two together using violin plots.
bind_rows(
"Real topics" = new_data,
"Fake topics" = new_data_fake %>%
rename(
topic_part = chunk_part,
speaker_topic = speaker_chunk
),
.id = "source"
) %>%
mutate(
# Make sure chunks are in correct order in plot.
topic_part = factor(topic_part, levels = c('start', 'middle', 'end'))
) %>%
ggplot(
aes(
x = topic_part,
y = prediction,
colour = source
)
) +
geom_violin(
draw_quantiles = c(0.25, 0.5, 0.75),
alpha = 0.2
) +
labs(
title = "Predicted Values for Real and Fake Topics at Midpoint of Monologue",
y = "Predicted (scaled) amplitude",
x = "Part"
)
Figure 3.17: Predicted values of real and fake topics at mid point of monologue.
While the violin plots in Figure 3.17 might not look very different from one another, we see the increase over the mean amplitude at the start of the real topics and a greater drop at the end of the monologue. The fake topics, by contrast, stay much closer to the mean value for amplitude. The reduction from start to end in the real topics is great than the reduction in the fake topics.
We can make this more rigourous by estimating standard errors for both models
using the confint.merMod function. This takes quite a while to compute. So
we save the results and load them in the next cell. To run this yourself,
change eval to TRUE, as above.
boot_real <- confint.merMod(
glmm_fit,
method = "profile",
parm=c(
'(Intercept)',
'topic_partmiddle',
'topic_partend',
'speaker_scaled_time',
'speaker_scaled_pitch'
)
)
boot_fake <- confint.merMod(
glmm_fit_2,
method='profile',
parm=c(
'(Intercept)',
'chunk_partmiddle',
'chunk_partend',
'speaker_scaled_time',
'speaker_scaled_pitch'
)
)
write_rds(boot_real, here('models', 'real_topic_boot.rds'))
write_rds(boot_fake, here('models', 'fake_topic_boot.rds'))
boot_real <- read_rds(here('models', 'real_topic_boot.rds'))
boot_fake <- read_rds(here('models', 'fake_topic_boot.rds'))
print('Real topics')
## [1] "Real topics"
boot_real
## 2.5 % 97.5 %
## (Intercept) 0.03687131 0.08162179
## topic_partmiddle -0.03536956 0.01507101
## topic_partend -0.10206502 -0.04170599
## speaker_scaled_time -0.29140439 -0.11502071
## speaker_scaled_pitch 0.53472079 0.54867565
print('Fake topics')
## [1] "Fake topics"
boot_fake
## 2.5 % 97.5 %
## (Intercept) 0.04564677 0.09253894
## chunk_partmiddle -0.02397603 0.02818460
## chunk_partend -0.04271154 0.01917069
## speaker_scaled_time -0.26551632 -0.05152494
## speaker_scaled_pitch 0.54219642 0.55653181
Let’s plot these:
real_coefs <- coef(summary(glmm_fit))[, 'Estimate'] %>%
as_tibble(rownames = "variable") %>%
mutate(
ll = boot_real[,'2.5 %'],
ul = boot_real[, '97.5 %']
)
fake_coefs <- coef(summary(glmm_fit_2))[, 'Estimate'] %>%
as_tibble(rownames = "variable") %>%
mutate(
ll = boot_fake[,'2.5 %'],
ul = boot_fake[, '97.5 %'],
variable = str_replace(variable, 'chunk', 'topic')
)
glmm_coefs <- bind_rows(
"Real" = real_coefs,
"Fake" = fake_coefs,
.id = "Source"
)
pd = position_dodge(width=0.5)
glmm_coefs %>%
filter(
!variable %in% c('speaker_scaled_time', 'speaker_scaled_pitch')
) %>%
mutate(
variable = factor(
variable,
levels = c(
'(Intercept)', 'topic_partmiddle', 'topic_partend'
)
)
) %>%
ggplot(
aes(
x = variable,
y = value,
colour = Source,
group = Source
)
) +
geom_line(position=pd, alpha = 0.5) +
geom_point(position=pd) +
geom_errorbar(
aes(
ymin = ll, ymax = ul
),
width = 0.25,
position=pd,
alpha = 0.5
) +
labs(
title = "Bootstrap 95% Confidence Intervals of Coefficients for Real and Fake Topics",
x = "Topic part",
y = "Scaled amplitude"
)
Figure 3.18: Bootstrap coefficients for real and fake topics.
The middle and end of the fake topics have 0 within their 95% confidence intervals. The reduction at the end in the real topics is outside the range to be expected from the fake topics.
We now preform the same analysis with GAMM models.
qb_gamm_fake_data <- chunks_filtered %>%
mutate(
speaker_chunk = str_c(Speaker, '_', chunk),
speaker_chunk = as.factor(speaker_chunk),
Speaker = as.factor(Speaker)
)
gam_fake_fit <- bam(
speaker_scaled_amp_max ~
s(chunk_time_scaled) +
s(speaker_scaled_time) +
s(speaker_scaled_pitch) +
s(chunk_time_scaled, speaker_chunk, bs = "fs", k=5, m=1) +
s(speaker_scaled_time, Speaker, bs = 'fs', k=5, m=1),
data = qb_gamm_fake_data,
method = 'fREML',
discrete = TRUE,
nthreads = 8
)
write_rds(gam_fake_fit, here('models', 'gam_fake_fit.rds'))
gam_fake_fit_summary <- summary(gam_fake_fit)
write_rds(gam_fake_fit_summary, here('models', 'gam_fake_fit_summary.rds'))
gam_fake_fit <- read_rds(here('models', 'gam_fake_fit.rds'))
gam_fake_fit_summary <- read_rds(here('models', 'gam_fake_fit_summary.rds'))
gam_fake_fit_summary
##
## Family: gaussian
## Link function: identity
##
## Formula:
## speaker_scaled_amp_max ~ s(chunk_time_scaled) + s(speaker_scaled_time) +
## s(chunk_time_scaled, speaker_chunk, bs = "fs", k = 5, m = 1) +
## s(speaker_scaled_time, Speaker, bs = "fs", k = 5, m = 1)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.017769 0.006925 2.566 0.0103 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(chunk_time_scaled) 1.003 1.005 0.249 0.617
## s(speaker_scaled_time) 5.109 5.968 5.094 3.42e-05 ***
## s(chunk_time_scaled,speaker_chunk) 855.518 2634.000 0.724 < 2e-16 ***
## s(speaker_scaled_time,Speaker) 356.935 764.000 3.253 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0858 Deviance explained = 10.1%
## fREML = 92146 Scale est. = 0.72372 n = 72351
Let’s have a look at the main effect smooths:
plot_smooth(gam_fake_fit, view="chunk_time_scaled")
## Summary:
## * chunk_time_scaled : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
## * speaker_scaled_time : numeric predictor; set to the value(s): 0.474838460403316.
## * speaker_chunk : factor; set to the value(s): QB_NZ_F_408_2. (Might be canceled as random effect, check below.)
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(chunk_time_scaled,speaker_chunk),s(speaker_scaled_time,Speaker)
##
Figure 3.19: Smooth over fake topics
The effect here is less clear than that in Figure 3.12. We could, for instance, draw a straight line through these confidence intervals at around a scaled amplitude of 0.055. Note that we are asking the model to assume that we are roughly half way through the monologue, so we should not be surprised to be slightly above 0 throughout ( see Figure 3.20).
plot_smooth(gam_fake_fit, view="speaker_scaled_time")
## Summary:
## * chunk_time_scaled : numeric predictor; set to the value(s): 0.500786516853932.
## * speaker_scaled_time : numeric predictor; with 30 values ranging from 0.000000 to 0.999441.
## * speaker_chunk : factor; set to the value(s): QB_NZ_F_408_2. (Might be canceled as random effect, check below.)
## * Speaker : factor; set to the value(s): QB_NZ_M_626. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(chunk_time_scaled,speaker_chunk),s(speaker_scaled_time,Speaker)
##
Figure 3.20: Smooth over speaker scaled time.
This is, in general, quite similar to (Figure ??). This is good, as we have not removed the temporal order of the data. All we have done is select random places at which to start and end topics.
Further experimentation with these models reveals that the process of randomly
assigning ‘chunks’ has a big effect on the results of the models here. One
way to get a sense of things which does not depend on the specific ‘chunks’ is
to rerun the analysis multiple times as we did in our permutation tests in
corpus_pca.Rmd.
The script which performs this analysis is available at
scripts/permutation_topics.R. It reruns the chunking process and then fits
the linear mixed model structure above. It then extracts the coefficients,
t values, and errors from each model. 1000 iterations of the process were
carried out.
We load in the results here:
coefs <- read_rds(here('processed_data', 'glmm_pitch_coeffs_perm.rds'))
ses <- read_rds(here('processed_data', 'glmm_pitch_ses_perm.rds'))
t_vals <- read_rds(here('processed_data', 'glmm_pitch_tvals_perm.rds'))
We will visualise the distribution of coefficient results, errors, and t values.
We first define a general function:
violin_three_way <- function(plot_data, comparison_data = list()) {
out_plot <- plot_data %>%
as_tibble() %>%
pivot_longer(
cols = everything(),
values_to = "value"
) %>%
mutate(
name = factor(
name,
levels = c(
'(Intercept)', 'chunk_partmiddle', 'chunk_partend',
'speaker_scaled_time', 'speaker_scaled_pitch'
)
)
) %>%
filter(
!name %in% c('speaker_scaled_pitch')
) %>%
ggplot(
aes(
x = name,
y = value
)
) +
geom_violin(draw_quantiles = c(0.05, 0.5, 0.95))
if (length(comparison_data) > 1) {
out_plot <- out_plot +
geom_point(colour = "red", data = comparison_data)
}
out_plot
}
We then extract the coefficients, t-values, and ses from our model fit on the real data.
var_names <- c(
'(Intercept)', 'chunk_partmiddle', 'chunk_partend', 'speaker_scaled_time',
'speaker_scaled_pitch'
)
summ_glmm_fit <- summary(glmm_fit)
real_model <- summ_glmm_fit$coefficients[,1] %>%
as.tibble() %>%
mutate(
name = var_names
) %>%
filter(
!name %in% c('speaker_scaled_pitch')
)
real_model_s <- summ_glmm_fit$coefficients[,2] %>%
as.tibble() %>%
mutate(
name = var_names
) %>%
filter(
!name %in% c('speaker_scaled_pitch')
)
real_model_t <- summ_glmm_fit$coefficients[,3] %>%
as.tibble() %>%
mutate(
name = var_names
) %>%
filter(
!name %in% c('speaker_scaled_pitch')
)
We visualise the coefficients:
violin_three_way(coefs, real_model) +
labs(
title = "Distribution of Coefficient Estimates for Random and Topical Segments",
caption = "Red points indicate values for topical segments",
y = "Speaker scaled max amplitude",
x = "Coefficient"
)
Figure 3.21: Distribution of Coefficient Estimates for Random and Topical Segments.
The coefficient estimates for our real model are easily within the range
estimated in the ‘chunks’. Interestingly, on the whole, estimates for
chunk_partmiddle are lower than 0, and lower than the initial and final
segments of chunks. This the average model of our faked chunks will have a
‘u’ shape rather than a steady decline.
violin_three_way(t_vals, comparison_data = real_model_t) +
labs(
title = "Distribution of t-values for Random and Topical Segments",
caption = "Red points indicate values for topical segments",
y = "t-value",
x = "Coefficient"
)
Figure 3.22: Distribution of t-values for Random and Topical Segments.
Figure 3.22 that our real data sits towards the tails of the significance
achieved by the models fit on fake topics. We should be particularly surprised
by the values for chunk_partend and (Intercept), that is the first third of
the chunk. We see that the vast majority of the models we fit will have t-values
less than 2 for their various coefficients.
violin_three_way(ses, comparison_data = real_model_s) +
labs(
title = "Distribution of Standard Errors for Random and Topical Segments",
caption = "Red points indicate values for topical segments",
y = "ses value",
x = "Coefficient"
)
Figure 3.23: Distribution of standard errors for random and topical segments.
In all cases, we see the standard errors for the real topics sitting below the distribution for the fake topics. There are two phenomena which might be behind this:
We take it that the above provides at least some evidence that amplitude indicates position within a topical subsection of a monologue. The difficulty in pinning this down with increased confidence is the difficulty of distinguishing between the reduction in amplitude across the monologue as a whole and the reduction in amplitude which might be explained by position in a topical subsection. Straight-forwardly, the end of a topical section is later in the monologue than the start and so would be expected to be of somewhat lower amplitude.
See preprocessing.Rmd.↩︎
Wikipedia has detailed accounts of the September, 2010 and February, 2011 earthquakes.↩︎
Given that this is exploratory research, surprisingness is all we can get out of statistical test values (see Baayen et al. 2017, 227).↩︎